百山祖冷杉叶片转录组分析

王光炯，柳新红，许大明，吴义松，李因刚，于明坚，张丽芳，*

（1.浙江省林业科学研究院，浙江杭州 310023；2.浙江农林大学林业与生物技术学院，浙江杭州 311300；3.浙江凤阳山百山祖国家级自然保护区管理局百山祖管理处，浙江庆元 323800；4.浙江大学生命与科学学院，浙江杭州 310058）

【研究意义】百山祖冷杉（Abies beshanzuensisM.H.Wu）为松科（Pinaceae）冷杉属（Abies）常绿乔木[1]。我国冷杉属植物种类丰富，约有22个种3个变种，其中如百山祖冷杉、元宝山冷杉（Abies yuanbaoshanensisY.J.Lu et L.K.Fu）、梵净山冷杉（Abies fanjingshanensisW.L.Huang，Y.L.Tu&S.Z.Fang）和资源冷杉（Abies beshanzuensisvar.ziyuanensis（L.K.Fu et S.L.Mo）L.K.Fu et Nan Li）等分布狭窄且数量极少，为国家一级重点保护植物[2]。百山祖冷杉于1963年由吴鸣翔先生首先发现，1976年正式命名发表，打破了许多学者认为我国东南大陆地区不可能有冷杉分布的论断[3-4]。它自然分布于浙江省丽水市庆元县百山祖国家级自然保护区海拔1 740～1 750 m处，目前原生地仅存3株野生成熟植株。现阶段由于百山祖冷杉现存数量稀少，难开花结果，自然繁殖十分困难等原因，急需加强科学监测、研究和有效保护[8-10]。在百山祖冷杉保护和研究中对濒危机制有初步的了解，除了开“花”失败、自身花粉无效传播和无法产生有活力种子等原因以外，幼苗和幼树由于环境的变化（气候变暖和环境污染）等原因不断地死亡。然而，到目前为止，对于百山祖冷杉的生物学-生态学特性认识严重不足，适应胁迫环境的分子调控机制也尚不清楚。【前人研究进展】伴随后基因组时代来临，转录组学连接着基因组与蛋白质组，是率先发展并应用最广、最活跃的一门技术。近年来，国内外的学者对许多濒危植物展开了转录组测序研究，例如，在建立多叶斑叶兰（Goodyera foliosa（Lindl.）Benth.）的繁殖体系的研究中同时结合高通量转录组测序技术，成功建立植株再生体系并得到全面完整的转录组信息，为其扩繁、基因功能、遗传发育以及调控机制研究奠定基础[5-7]。【本研究切入点】目前Illumina、Pacbio、SLAF-seq等技术被广泛应用于功能基因研究、分子标记的开发、代谢途径研究以及应对各种环境胁迫反应机制研究等方面[11-14]。在冷杉属植物朝鲜冷杉（Abies koreana）的研究中，通过二代转录组分析技术分析热胁迫下表达的基因，提供了对遭受热胁迫的朝鲜冷杉的第一个综合表征，为识别耐高温胁迫的品系提供了重要的参考指标[15]。【拟解决的关键问题】本文基于高通量二代转录组测序平台对百山祖冷杉的转录组进行测序，获得大量的基因数据，并对组装的转录本进行分类统计、功能注释，通过分析转录组数据研究百山祖冷杉代谢途径关键基因，为后续的环境胁迫调控机理等方面提供研究基础。另外，通过转录组测序和组装，开发微卫星分子标记，研究百山祖冷杉母树和子代个体的遗传多样性和亲缘关系。

1 材料与方法

1.1 试验材料

于2019年6月在浙江省庆元百山祖保护区内（27°45′N，119°11′E）采集百山祖冷杉测序材料，选取当年萌发枝条上部幼嫩叶片，3个生物学重复，采集后放置于干冰内短期保存，之后保存于-80 ℃超低温冰箱，用植物RNA提取试剂盒提取RNA。

1.2 试验方法

1.2.1 RNA提取、文库构建以及Illumina测序 从百山祖冷杉叶片组织中提取总RNA，并且检测RNA的纯度和浓度[16-17]。构建基于Illumina Hiseq高通量测序的cDNA文库。文库构建后，本研究为无参百山祖冷杉转录组分析，先拼接由测序所得到的序列，以拼接好的转录本为参考序列，利用GO、KO、KOG、NR、等7大数据库与转录组测序结果比对并进行功能注释，再分析其GO富集和KEEG富集[18]。

1.2.2 基因表达量分析 利用BioKanga‘maploci’将读长进行FPKM（Reads per kilobase per million reads）计数标准化[19]。FPKM能消除基因长度和测序量差异对计算基因表达的影响。利用Bowtie软件将Reads与Unigene库进行比对并结合RSEM软件对表达量水平进行评估，通过FPKM值来表示对应Unigene的表达丰度。

1.2.3 基因功能注释 通过BLAST软件将所得到的Unigene序列与NR、Swiss-Prot、GO、COG、KOG、egg‐NOG4.5、KEGG 7大公共数据库进行比对，得到所需的Unigene注释信息[20]。

1.2.4 转录因子分析 将测序所得到的基因序列与转录因子数据库在BLAST软件中进行比对[21]。

1.2.5 SSR位点筛选和引物设计 通过MISA（http:∕∕pgrc.ipkgatersleben.de∕misa∕）软件对筛选得到的1 kb以上的Unigene进行SSR分析，根据碱基的重复次数鉴定SSR类型[22]。通过Primer 3引物设计软件对找到的SSR位点序列进行引物设计[23]。通常引物设计需要注意的有以下几点:（1）扩增的产物大小在100～300 bp，所设计引物的长度尽量在18～25 bp，引物序列中的GC含量控制在40%～65%；（2）退火温度56～64 ℃，并且设计的上下游引物退火温度的差值小于5 ℃。设计引物时，引物不能出现SSR位点序列。（3）设计时应避免出现发卡结构、二聚体以及错配等情况。随机选取10对引物，委托上海捷瑞生物有限公司合成。

2 结果与分析

2.1 百山祖冷杉转录组测序结果以及数据组装

基于二代测序技术Illumina Hiseq高通量测序平台对百山祖冷杉的叶片进行转录组测序，最终获得大量原始数据（Raw data），原始测序序列去除一些低质量以及含有测序引物、接头等人工序列的reads后，共获得7.02 Gb Clean Data，Q30碱基百分比在95.43%及以上。利用Trinity软件对Clean Data进行组装，组装得到73 103条转录本序列，总长度为9 303 348 bp，N50长度为1 930 bp，平均长度为1 272.64 bp。得到41 228条Unigene，总长为42 458 236 bp，Unigene的N50为1 690 bp，平均长度1 029.84 bp，组装完整性较高（表1）。

表1 拼接长度分布Tab.1 Splicing length distribution

2.2 百山祖冷杉转录组基因功能注释

通过BLAST（E-value＜1e-5）软件将百山祖冷杉测序所得到的所有Unigene与Nr、eggNOG4.5、Swiss-Prot、KOG、COG、GO、KEGG数据库进行比对，得到有注释信息的Unigene为27 738个，占所有Unigene的67.28%。其中注释到Nr数据库的数量最多为27 418，占所有Unigene的66.50%。注释到COG数据库的Unigene数量最少，仅有9 411条占所有Unigene的22.82%（表2）。

表2 Unigene注释统计Tab.2 Unigene annotation statistics

Nr数据库中E值分布的分析结果显示，E值分布于0～1e-150的Unigene最多，高达6 711条，占注释到Nr数据库的Unigene的24.48%。同一性分布显示，匹配的同一性为80%～90%的Unigene最多，达到4 908条，占17.90%。另外，对Nr注释到基因同源序列的物种分布进行了分析（图1）。

2.3 Unigene的功能分类

将百山祖冷杉的全部Unigene与GO数据库进行比对，并将比对到的功能分类进行统计。根据结果表明有16 020条Unigene在GO数据库中注释，匹配包括了3大类，分别为细胞组分、分子功能和生物学过程。注释到细胞组分大类共15个亚类的基因有33 892个，注释到分子功能大类共15个亚类的数量为18 407个，注释到生物学过程大类共20个亚类的基因共有29 809个。在细胞组分的大类中，细胞（cell）和细胞组分（cell part）占比最高，分别为6 891和6 890个。在分子功能的大类中催化活性（catalytic activity）亚类注释的数量最多为8 185个，其次绑定功能（binding）注释数量为7 349。生物学过程大类中占比最高的亚类为代谢过程（metabolic process）注释数量为8 080，其次细胞过程（cellular process）注释数量为7 107（图2）。

图1 Unigene在Nr数据库中的各项分布Fig.1 Distribution of Unigene in Nr Database

图2 Unigene的GO功能分类Fig.2 GO function classification diagram of Unigene

将百山祖冷杉的所有Unigene与COG数据库进行比对，发现有9 411条Unigene得到注释。分配至26类功能中，基因的数量在不同的功能分类存在明显的差异。其中，仅一般功能预测类基因数量最多，达到1 261条，其次糖运输和新陈代谢功能类别，注释的基因多达1 028条。此外，翻译、核糖体结构、生物合成功能类别（970条）与翻译后修饰、蛋白质折叠和分子伴侣功能类别（813条）注释得到的基因也较多（图3）。

图3 COG功能分类统计Fig.3 COG function classification statistics

2.4 基因转录因子分析

在植物的生长发育过程中常会遇到干旱、低温、洪涝等多种非生物因子的胁迫作用。在植物的正常生长发育中常常会遭受逆境的影响，影响严重时甚至对植物的伤害是不可逆的，最终导致植株的死亡[24]。植物中含有大量的转录因子，当植物参与逆境反应时相关转录因子调控胁迫应答基因表达，使植物产生抗逆反应，通过生理生化的变化来减少自身在逆境中受到的伤害。百山祖冷杉转录组测序共获得1 872个转录因子，分布在208个转录因子家族中（图4）。其中C2H2、AP2∕ERF-ERF、MYB-related、bHLH、MYB、C3H等转录因子家族所包含的转录因子较多，分别占4.70%、3.63%、3.25%、2.99%、2.93%、2.88%，这些转录因子家族主要对植物生长发育的调控起到重要作用。

图4 转录因子家族分布Fig.4 Transcription factor family distribution

表3 部分转录家族及相应Unigene数量Tab.3 Partial transcription factor family and corresponding Unigene numbers

图5 Unigene的代谢通路分析Fig.5 Unigene’s metabolic pathway analysis

2.5 Unigene的代谢通路分析

将所得到的Unigene与KEGG数据库进行比对后进行代谢通路的分析，结果显示共有9 976条Unigene成功获得注释。根据功能将这些所注释的Unigene分为6个大类，其中有21个亚类涉及126条代谢途径。6大类分别为新陈代谢、遗传信息处理、细胞过程、有机系统、环境信息处理、人类疾病。新陈代谢大类涉及95条代谢通路3 102个基因，遗传信息处理大类涉及21条代谢途径1 949个基因，细胞过程大类涉及3条代谢途径471个基因，环境信息处理大类涉及3条代谢途径272个基因，有机系统和人类疾病大类均涉及2条代谢途径所包含的基因分别为162个和15个。其他代谢通路的详细信息见图5。

2.6 百山祖冷杉转录组SSR位点分析

将所得到的数据进行筛选，选出长度大于1 kb的Unigene，通过MISA软件进行SSR分析，共得到1 794个SSR位点，相比其他物种的SSR位点数量，百山祖冷杉的位点数量显著减少。利用Primer 3软件对SSR位点两侧的序列进行引物设计，共获得100多对SSR引物，部分引物信息见表4。

表4 百山祖冷杉SSR部分引物信息Tab.4 Information of part pairs of development from Abies beshanzuensis

2.7 百山祖冷杉基因表达量的分析

百山祖冷杉的转录组测序中共得到23 540 099条reads，其中有20 271 952条reads比对到基因图上，占总量的86.12%，有13 041 470条reads比对到基因图唯一位置，占总量的64.33%，比对到多个位置的reads为7 230 482占35.67%。测序得到73 103条转录本序列，41 228条Unigene。41 228条Unigene的FPKM值的平均值为17.28，范围在0～28 221.41，其中，6 786条Unigene的FPKM值在平均值以上。

3 结论与讨论

生物多样性丧失和物种灭绝是当今世界面临的重大生态危机之一，许多珍稀濒危植物是自然生态系统中的关键种，在维持自然生态系统运行中发挥着重要作用，一个物种的灭绝，有可能发生连锁反应，对自然生态系统的稳定性造成一定的影响，导致灾难性后果[25-27]。百山祖冷杉从1976年发现以来，基础研究取得了一些成果，但是在外界环境中，百山祖冷杉对环境因子的胁迫反应机制一直不是很清楚。在植物的逆境胁迫中，经常受到非生物逆境（如干旱、盐、极端温度及重金属胁迫等）的侵害，在应答干旱胁迫方面，DREB、WRKY、NAC、bHLH及bZIP等转录因子受干旱信号诱导，调节下游抗旱基因的表达，从而提高蔬菜作物抗旱能力。同时，水分运输相关功能基因（PIP、TIP）、E3连接酶SIZ1及脱水蛋白DHN也被报道受干旱诱导，并通过调节水势、渗透势及ROS积累抵御干旱胁迫。在抵御盐胁迫方面，SOS途径至关重要。Sl SOS2能够通过调节Sl SOS1和Na+∕H+逆向转运蛋白Le NHX2∕4的表达维持离子平衡和调节植物器官中Na+的分配。蔬菜抗盐研究中NAC、ERF、MYB等转录因子响应盐胁迫并激活抗逆相关基因表达，从而提高蔬菜作物抗盐能力[28]。

本文利用RNA-seq技术[29]，对百山祖冷杉的转录组进行测序、组装，并与Nr、GO、KEGG等7大公共数据库进行比对。GO分析显示，获得注释的基因共参与了43个通路，最富集的通路主要是在生物学过程中，包括免疫系统形成、发育过程和代谢等，而代谢过程获得的注释是最多的。而在COG数据库比对分析结果中，糖运输和新陈代谢功能注释的基因最多，植物在干旱胁迫条件下，体内水分代谢与碳代谢会发生失衡现象，光合速率降低、蒸腾速率降低，导致生长降低，使植物陷入饥饿状态[30]。糖作为重要的代谢物质，生物合成不仅为植物的生长发育提供碳源和能源，而且分解产生的小分子单糖可以调控植物细胞的渗透压以抵抗胁迫[31]。GO分析和COG数据库分析的结果一致，进一步说明百山祖冷杉在生长过程中代谢活动旺盛，而且一些次生代谢物的生物合成途径，包括生物碱、萜类化合物和黄酮类化合物等在抵抗病原物的入侵方面起着重要的作用[32]。

转录因子在植物的生长和调控中起着重要的作用[33]。在植物对抗环境胁迫过程中，调控相关生理反应基因表达的蛋白在防卫反应和逆境信号转导中发挥作用，使植物适应外界不良环境[34-36]。百山祖冷杉转录因子分析表明，C2H2、AP2∕ERF-ERF、MYB-related、bHLH、MYB、C3H等转录因子家族所包含的转录因子较多，这些转录因子家族在植物发育中起到关键调控作用。已有的研究表明，C2H2锌指转录因子家族主要参与植物叶的发生，花器官的调控，侧枝的形成及逆境胁迫等生命过程[37]，意味着C2H2锌指转录因子家族在百山祖冷杉生长各阶段、各组织中都参与表达。另外，有研究表明AP2∕ERF-ERF家族成员能够参与调控植物花的形成和发育，种子的大小，种子的重量，种子中脂肪酸的合成和抗逆等。该家族基因在种子不同部位存在特异性表达，在种子发育的各个过程中起着不同的作用[38]。各转录因子分析研究为百山祖冷杉种子败育的机理研究提供参考，同时为今后的幼苗繁育及生境改良工作提供分析参考。

通过转录组测序得到的数据开发SSR引物，应用于百山祖冷杉遗传多样性等方面研究。前人研究中已有对百山祖冷杉遗传多样性的研究[39-40]，但大多为ISSR、RAPD等分子标记[38]，遗传多样性研究中存在引物数量有限等原因，无法将百山祖冷杉的杂交实生苗或无性繁殖的个体与其亲本进行鉴别[41]。通过转录组的高通量测序对数据进行筛选和SSR分析得到1 794个SSR位点，所得到的位点数量并不多，重复类型以单核苷酸重复、双核苷酸重复和三核苷酸重复为主，其中单核苷酸重复类型最多。根据这些位点信息设计了大量的SSR引物，为今后百山祖冷杉母本与子代亲缘关系鉴别以及进一步的遗传多样性分析提供基础信息。

本研究基于百山祖冷杉转录组测序所获得的大量数据，通过后续分析系统了解了百山祖冷杉的基因功能注释、代谢通路等，为今后对其逆境胁迫、遗传背景、基因功能以及遗传标记等研究提供基础数据，为百山祖冷杉的保护以及种群的扩大繁殖提供理论基础。但现阶段对于百山祖冷杉的研究刚在起步阶段，由于技术手段的限制，未对研究中所注释的Unigene基因进行定量PCR验证，进一步确定转录组数据的准确和可实用性，而这些存在的问题也是接下来百山祖冷杉研究工作的重点。