濒危植物秦岭石蝴蝶花瓣数量变异机理研究

蒋景龙，孙 旺，李 丽，李 耘，胡凤成

(1 陕西理工大学 生物科学与工程学院，陕西汉中 723000；2陕西理工大学 化学与环境科学学院，陕西汉中 723000； 3 汉中市野生动植物保护管理站，陕西汉中 723000；4 略阳县苗圃，陕西汉中 724310)

花是被子植物特有的观赏和生殖器官，是植物分类学中的重要依据，同时也是被子植物进化过程中变化最丰富的器官[1]。研究花器官发育及其分子机制对其植株的发育及演化过程具有重要意义[2-3]。花器官发育是由相关基因经过一系列复杂的调控机制形成的，Coen等[4]首次于1991年提出了花器官发育的经典ABC模型以及随着研究的深入，由ABC模型发展起来的其他模型。基于模式植物拟南芥和金鱼草各类花器官的研究，在其他物种中大量决定花器官特征属性的基因相继被克隆鉴定出来，任何一种基因突变都会导致花器官形态发生变化[5-6]。在环境因素(不同生物和非生物因素)的影响下，花器官形态常会显现出特定的变异类型，其分子机制归于基因结构的改变和表达模式的多样化[7-8]。对花器官变异的研究有助于揭示濒危种群野生植物响应环境变化的遗传机制，对理解其致濒机理进而采取相应的保护措施具有重要的意义[9]。

秦岭石蝴蝶(PetrocosmeaqinlingensisW. T. Wang)属苦苣苔科(Gesneriaceae)石蝴蝶属(PetrocosmeaOliv.)中华石蝴蝶组(Sect.Petrocosmea)的多年生草本植物，其淡紫色的花冠和翠绿覆绒毛的叶片具有很高的观赏价值。目前该物种仅在陕西省汉中市勉县和略阳县境内的阴湿山沟发现，分布区域狭窄，数量稀少，且极易受环境和人为活动的影响，被列为秦岭地区特有的国家二级珍稀濒危保护植物[10-12]。本课题组通过无性繁殖技术人工繁殖秦岭石蝴蝶10 000余株，在人工繁育的这些秦岭石蝴蝶群体中发现了大量的花器官变异，除中国植物志中描述的正常秦岭石蝴蝶花冠为二唇形，两侧对称，上唇2裂、下唇3裂形态外，共发现了17种变异的花冠类型[13]。

本研究采用Illumina HiSeq 2500高通量测序技术对秦岭石蝴蝶正常花朵2-3型(上唇2裂、下唇3裂，占65.43%)和在所有的变异类型中比例最多的变异类型2-4型(上唇2裂、下唇4裂，占26.15%)的开花早期和晚期进行了转录组测序，为解析秦岭石蝴蝶的花器变异调控机制提供一定的参考。

1 材料和方法

1.1 植物材料采集及形态学观察

实验材料于2019年7月2日采自略阳县苗圃的人工仿野生栽培大棚中，所取秦岭石蝴蝶均来自同一苗床，株龄、植株大小、生长状态均相对一致。选取处于发育阶段正常花朵类型(2-3型，N)和典型变异花朵类型(2-4型，V)的花苞(A)和成熟开放的花朵(B)为实验材料。实验材料分为4组：正常花苞(NA)、变异花苞(VA)、正常花朵(NB)和变异花朵(VB)，每组3个重复分别编号。保存在干冰中带回实验室，-80 ℃保存，选取完整的正常花和变异花花苞和花朵，拍照保存。

1.2 方 法

1.2.1 秦岭石蝴蝶花的总RNA提取、文库构建及转录组测序秦岭石蝴蝶花总RNA提取方法参照Trizol Reagent方法进行，用RNA专用琼脂糖电泳检测RNA的浓度和纯度。通过Oligo(dT)磁珠富集总RNA中带有polyA结构的mRNA，采用离子打断的方式，将mRNA打断到300 bp。以片段mRNA为模板，用6碱基随机引物和逆转录酶合成cDNA第一链，并以第一链cDNA为模板合成第二链cDNA，文库构建完成后，采用PCR扩增进行文库片段富集，之后根据片段大小进行文库选择，文库大小在450 bp。通过Agilent 2100 Bioanalyzer对文库进行质检，再对文库总浓度及文库有效浓度进行检测，形成单链文库。建库后，由上海派森诺公司完成测序，测序平台为Illumina HiSeq。

1.2.2 数据处理、转录本拼接及功能注释首先对原始下机的FASTQ文件数据进行过滤，将带接头、长度小于50 bp、序列平均质量在Q20以下的Reads进行去除，使用Trinity软件对高质量序列进行拼接得到转录本序列(Transcript)，挑选每个基因下最长的转录本作为该基因的代表序列Unigene。对Unigene与六大数据库进行基因功能注释，基因功能注释所用到的数据库包括NR(NCBI non-redundant protein sequences)、GO(Gene Ontology)、KEGG(Kyoto Encyclopedia of Genes and Genome)、eggNOG(evolutionary genealogy of genes: Non-supervised Orthologous Groups)、Swiss-Prot和Pfam。

1.2.3 表达分析使用R语言的DESeq软件包，根据表达量对各样品进行PCA主成分分析。采用DESeq对基因表达进行差异分析，筛选差异表达基因条件为：表达差异倍数|log2FoldChange|>1，显著性P<0.05。对不同分组之间的差异表达基因进行统计，根据差异分析的结果，统计各比较组之间的共有和特有差异基因数量。

1.2.4 差异基因富集分析使用topGO进行GO富集分析，分析时利用GO term注释的差异基因对每个term的基因列表和基因数目进行计算，然后通过超几何分布方法计算P(显著富集的标准为P<0.05)，找出与整个基因组背景相比，差异基因显著富集的GO term，从而确定差异基因行使的主要生物学功能。Pathway注释主要采用KEGG的KAAS(KEGG Automatic Annotation Server)在线自动化注释系统完成(http://www.genome.jp/tools/kaas/)。

1.2.5 差异表达基因的实时荧光定量PCR分析为了验证基因差异表达的测序结果，选取差异表达的基因进行荧光定量PCR验证。用SPSS19.0软件进行生物统计学分析，通过双侧t检验方法进行分析，P<0.05为差异有统计学意义。

2 结果与分析

2.1 两个时期正常(2-3型)和典型变异(2-4型)花型分析

在解剖秦岭石蝴蝶正常小花时发现，花器官形成早期上唇和下唇出现明显的深裂，形成了早期的上唇2裂，下唇3裂(图1，A)，成熟花朵花冠二唇形，上唇2裂，下唇3裂，命名为2-3型(图1，B)；变异小花早期出现了上唇2裂，下唇4裂的深裂(图1，C)，最终形成了变异类型中的上唇2裂，下唇4裂，简称2-4型(图1，D)。

2.2 秦岭石蝴蝶的花瓣数量变异的转录组分析

2.2.1 样本基础数据分析对12个样本上机测序得到的原始数据分别进行筛查和过滤，去除带接头、低质量的序列，获得高质量序列34.55～47.27 Mb个，并统计其占原始序列的比例，占比为88.44%～93.42%，高质量序列碱基数为5.18～7.09 Gb，GC含量为42.78%～44.56%(表1)。对高质量序列进行拼接共获得Unigene序列113 889个，其总长度为99 295 924 bp，平均长度872 bp，N50长度为1 592 bp，GC含量为38.57%(表2)。

A.正常2-3型花苞(NA)；B.正常2-3型花朵(NB)；C.变异2-4型花苞(VA)；D.变异2-4型花朵(VB); a.上唇瓣；b.下唇瓣；c.花药；d.子房；NA、NB、VA、VB下同图1 秦岭石蝴蝶正常2-3型和变异2-4型花苞和花朵A. Normal 2-3 type flower bud(NA); B. Normal 2-3 type flower(NB); C. Variation 2-4 type flower bud(VA); D. Variation 2-4 type flower(VB); a. Upper lip; b. Lower lip; c. Anther; d. Ovary； NA, NB, VA and VB are same as belowFig.1 Flower buds and flowers of normal type 2-3 and variation type 2-4 of Petrocosmea qinlingensis

2.2.2 秦岭石蝴蝶花器官测序比对本研究将秦岭石蝴蝶的Unigene与六大功能数据库进行比对，并根据基因的相似性进行基因功能注释，最终分别有52 677(NR: 46.25%)、46 750(eggNOG: 41.05%)、35 069(Swissprot: 30.79%)、21 088(GO: 18.52%)、19 663(KEGG: 17.27%)和18 215(Pfam: 15.99%)个Unigene获得功能注释，其中在所有数据库中都被注释到的Unigene数目有5 367个。通过与NR库的同源序列相比，E值较低(1×10-30≤E<1×10-5)的Unigenes占Unigenes总数的42.81%，E值极低(1×10-100≤E<1×10-30)的Unigenes占Unigenes总数的35.19%，表明注释的Unigenes序列与NR数据库的同源序列具有高度相似性。

与NR库进行比对注释，可以获取本物种基因序列与近缘物种基因序列的相似性信息。在NR数据库中，秦岭石蝴蝶与旋蒴苣苔(Dorcocerashygrometricum) Unigene序列匹配度最高，占注释Unigene总数的54.29%，其次是芝麻(Sesamumindicum)，占比10.38%，最低的是甜橙(Citrussinensis)，占比0.89%，还有22.21% Unigene与其他物种序列无匹配(图2)。

2.2.3 主成分分析与基因差异表达分析对12个样本进行PCA主成分分析(图3)，NB与VB距离较近，且组内分布较为集中，而NA与VA虽然也可分为一组，但组内分布较为分散。总体而言，花苞与花朵这两个不同发育时期基因表达差异明显，但4组材料基本分开，可以进行后续分析。采用DESeq对基因表达进行差异分析，统计各比较组之间的差异基因数量。韦恩图(图4)展示了各比较组之间的共有和特有差异基因数量。其中4个比较组中所共同拥有的差异基因数为62个，NA vs NB、VA vs VB、NA vs VA与NB vs VB独有的差异基因数量分别为3 079、2 914、676与1 382个。NA vs NB与VA vs VB差异表达基因分别为10 665(5 248个上调表达基因和5 417个下调表达基因)和10 608个(5 075个上调表达基因和5 533个下调表达基因)，明显高于NA vs VA与NB vs VB的2 917个(1 829个上调表达基因和1 088个下调表达基因)与3 500个(1 721个上调表达基因和1 779个下调表达基因)(表3)。

表1 样本测序数据统计

表2 转录组测序的数据序列统计

图2 NR数据库E值分布图和物种分布图Fig.2 E-value distribution and species distribution of NR Database

图3 12个样本的PCA分析Fig.3 PCA analysis of 12 samples

2.3 秦岭石蝴蝶花瓣数量变异的差异代谢通路分析

对差异表达基因进行GO富集分析，按照分子功能、生物过程和细胞组分进行GO分类，对每个比较组GO分类中p-value最小即富集最显著的前20个GO条目，取-lg(p-value)进行展示(表4)。其中NA vs VA差异基因在每个GO分类中富集最显著的GO条目分别为细胞膜、跨膜转运和反向转运蛋白活性，NB vs VB差异基因在每个GO分类中富集最显著的GO条目分别为酶抑制剂活性、细胞膜和催化活性的负调控，NA vs NB和VA vs VB差异基因在每个GO分类中富集最显著的GO条目均为细胞膜、细胞壁组织或生物合成和催化活性。

图4 组间差异表达基因韦恩图Fig.4 Venn diagram of differentially expressed genes between groups

对差异表达基因进行KEGG通路富集分析，将排名前20的KEGG pathway进行展示(图5)。通过富集因子(Rich factor)、FDR值和富集到此pathway上的基因个数来衡量富集的程度。Rich factor越大，表示富集的程度越大。FDR一般取值范围为0～1，越接近于零，表示富集越显著。NA vs VA差异基因主要富集在植物激素信号转导和脂肪酸延长等途径中；NB vs VB差异基因主要富集在戊糖和葡萄糖醛酸的相互转化、苯丙烷类生物合成、硫代谢、类黄酮生物合成、玉米素的生物合成等途径中；NA vs NB差异基因富集最显著的KEGG pathway主要有植物激素信号转导、植物-病原菌互作、氧化磷酸化、苯丙烷类生物合成、戊糖和葡萄糖醛酸的相互转化及淀粉和蔗糖代谢等；VA vs VB差异基因富集最显著的KEGG pathway主要有植物激素信号转导、苯丙烷类生物合成、戊糖和葡萄糖醛酸酯的相互转化及氧化磷酸化等。

表3 表达差异分析结果统计表

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

续表4 Continued Table 4

2.4 花瓣数量变异差异基因的表达分析

由于秦岭石蝴蝶目前无参考基因组，因此注释得到的许多基因功能未知或未预测，明确功能的基因很少。以校正后的P值为参数，将4个比较组(NA vs NB、VA vs VB、NA vs VA与NB vs VB)中的表达差异基因进行交叉比较，结合功能注释与富集结果，选出部分显著性强的表达差异基因。进一步对4个比较组进行交叉比对分析，获得6个可能与秦岭石蝴蝶花器官发育相关且功能已知的基因PqMIF2、PqMYB340、PqMYB305、PqGATA12、PqCCD4和PqZBED。

将筛选出的6个基因利用qRT-PCR对转录组测序结果进行验证(图6)，转录组结果中上调表达的基因其验证结果也为上调表达，转录组结果中下调表达基因其验证结果也为下调表达，表明基于转录组测序数据的基因差异表达分析结果是可信的。6个基因qRT-PCR分析得到的相对表达量中PqMIF2、PqMYB340、PqMYB305与PqCCD4在开放的花朵中表达量极为丰富，但在花苞组织中表达量极低；而PqGATA12则明显在花苞中表达量高于开放花朵中。此外，PqZBED在正常花器官材料中的表达量高于变异花器官。这与转录组表达谱分析趋势一致，说明基因表达谱的分析结果可靠。

图6 筛选出的6个基因在4个样品中的表达Fig.6 Expression of six genes screened in four samples

3 讨 论

花器官数目的变化可通过原基分裂、原基融合或次生原基的形成而产生，花原基细胞数目的增减可能导致花器官数目增减，而引起花器官变异却是基因和环境条件共同作用的结果[14]。在获得许可的情况下于2019年8月中旬对略阳县境内的野生秦岭石蝴蝶进行了野外观察和记录。遗憾的是在所有观察到的花朵中未发现变异现象，所有花朵与中国植物志中的记录相一致，而将采集到的野生秦岭石蝴蝶植株带回人工控制的环境和营养土进行栽培，发现部分植株在移栽后开花时开始表现出花瓣数量增多的现象。这表明，丰富的营养物质可能是导致秦岭石蝴蝶花器变异的主要原因。为了探明秦岭石蝴蝶花瓣数量变异原因，本研究采用Illumina HiSeq 2500高通量测序技术对秦岭石蝴蝶两种花型发育的早期和晚期进行转录组测序，挖掘参与其花发育相关的差异基因并探讨花器变异的可能机制。

高通量测序结果显示，秦岭石蝴蝶基因注释与旋蒴苣苔植物的基因数据库匹配度最高，达到了54.29%。秦岭石蝴蝶属于石蝴蝶属植物，旋蒴苣苔属于旋蒴苣苔属，二者同属于苦苣苔科，表明秦岭石蝴蝶的大部分基因测序结果可靠，可以用于后续基因功能分析。PCA和组间差异表达基因分析结果显示，秦岭石蝴蝶花器官发育早期和晚期的差异表达基因数量远远大于同一时期的正常花器官和变异花器官的差异表达基因数量。这表明，相对于整个花器官发育的过程中差异表达的基因，影响花器官向着某一方向变异的基因数目并不是很大。

GO富集分析结果显示，早期的花器官变异可能与蛋白质的跨膜转运、定位和细胞分裂素的合成等生物过程有关，例如转运蛋白、半乳糖-棉子半乳糖转移酶、Rho鸟嘌呤核苷酸交换因子等具有很高的活性，而后期花器官的变异可能与果胶、糖类、蛋白质的水解有关，例如内肽酶、肽酶抑制剂活性、核苷酸酶等有关，同时在这一过程中水杨酸的调控也起到了重要作用。无论是正常还是变异的花器官，他们在不同发育的时期都与细胞壁的形成和组装、糖代谢和氧化还原生物过程有关。而KEGG富集分析也进一步证实了这些结果。在花器官发育的早期变异过程中，植物激素信号相关的基因差异表达非常显著，而在后期的花器官变异过程中则表现出戊糖和葡萄糖类和苯丙烷类的代谢调节相关基因发生了差异表达。

值得关注的是植物激素信号转导通路差异基因富集最显著，推测这些差异基因可能参与多种激素代谢调控植物器官发育。深入研究关键基因所参与这些代谢的反应，对揭示秦岭石蝴蝶花器官变异的分子机制具有重要意义。MIF在花期结束过程中抑制细胞增殖，在花的发育过程中控制心皮的数量，MIF的表达不足可阻止配子体组织的形成[15]。此外，研究发现，拟南芥的Mini Zinc Finger (AtMIF2)蛋白和其在番茄(Solanumlycopersicum)中的同源蛋白SlIMA作为结合蛋白参与了花分生组织的终止，在控制心皮数量的发育阶段起到重要作用[16]。在本研究中MIF在秦岭石蝴蝶的花器官变异中表达比较低，但在花器官发育的调控过程中却表达很高，表明在引起花器官变异过程中此基因并未发挥重要的作用，但在整个花器官发育过程中发挥着重要作用。同样，转录因子MYB的变化趋势与MIF基因一致。GATA转录调节因子广泛存在于真核生物中。在调节植物生长、开花时间、花的发育以及光周期等过程中发挥着重要作用[17]。本研究中发现GATA基因无论在花器官早期还是后期的变异过程中均有非常高的表达量，可能是引起秦岭石蝴蝶花器官变异的一个重要基因，在后期的基因功能验证中值得关注。大量研究表明，类胡萝卜素裂解双加氧酶4(CCD4)在调节植物花色和果实颜色中具有重要作用[18]。在秦岭石蝴蝶转录组分析结果中，CCD4在开放后的花朵中过量表达，而在花苞中表达量较低，这与两者的颜色差异相一致。

对于秦岭石蝴蝶而言，这种人工条件下的花器官变异是有利的，增加的花器官数目可以增加其观赏性，在今后观赏价值的开发方面具有积极意义。本研究结果表明，环境改变可能在花器官发育的早期引起了某些花器官发育调控基因的改变，从而导致花器官发育和表达调控异常，但其分子机制仍需进一步研究、分析和验证。