长瓣兜兰花2个不同时期转录组分析

颜凤霞， 王莲辉， 田 凡， 李 涛

(1.贵州省林业科学研究院， 贵阳 550005； 2.贵阳市林草资源监测中心， 贵阳 550003)

长瓣兜兰(Paphiopedilumdianthum)为兰科兜兰属植物，分布于贵州、云南和广西等省区，因其独特的外观和优雅的花型而具有较高的观赏价值[1]。长瓣兜兰还是中国少有的多花兜兰，常受育种工作者的青睐，将其作为亲本进行育种，但由于遭受人类毁灭性采挖，野生资源极少，已处于濒危状态[2]。

王莲辉等[1]对长瓣兜兰的研究仅限于其栽培、组织培养及无菌播种快速繁殖体系等研究，在分子生物学层面的研究未见报道。兜兰属植物花发育及开花的分子机制及相关基因的挖掘也鲜见报道。RNA-seq 技术是近年来兴起用于在分子层面研究动物、植物、微生物的新技术[3-6]，该技术已应用于兰科植物多方面的研究。如在拟兰属中，Lin等[7]通过对其进行基因组测序，研究了其进化关系。Li等[8]通过转录组测序对扇叶文心兰中MADS-box家族基因进行分析。近年来已有两个兜兰品种通过转录组测序，开发出了SSR分子标记[9]。为了丰富兜兰分子生物学研究数据，挖掘与兜兰属花发育的相关基因。本研究利用RNA-seq 技术，对长瓣兜兰两个不同发育时期的花进行了转录组测序和分析，以期为解析长瓣兜兰花发育及开花的调控机制、培育优良兜兰品种奠定基础，并为其基因组测序组装提供参考。

1 材料与方法

1.1 实验材料

本研究以种植于贵州省林业科学研究院兰花资源圃中生长健壮的长瓣兜兰为实验材料， 7月和8月采集其花蕾和盛开花朵，每个样品三个生物学重复。花蕾和花朵采集后用液氮处理后于-80 ℃冰箱保存备用。

1.2 RNA提取与转录组测序

本研究用RNAprep pure植物总RNA提取试剂盒(天根)分别对长瓣兜兰花蕾和花朵的总RNA进行提取。然后用1%琼脂糖凝胶电泳检测总RNA的完整性，并用赛默飞公司的Multiscan GO对总RNA的浓度和纯度进行检测。

通过RNA 质量检测符合标准后，委托上海生工生物工程股份有限公司对长瓣兜兰花器官两个不同时期花进行转录组测序。

1.3 测序数据分析

测序得到的原始数据进行过滤，去除只含有测序接头序列和N 含量过高以及低质碱基(Q值<20)及过短的序列数据。利用Trinity软件拼接完成后，将所得核苷酸序列分别与NR、KEGG、KOG、Swissprot等数据库进行比对，获得与长瓣兜兰花器官转录本所对应的注释信息。使用MISA软件对基因转录本进行SSR检测。利用FPKM(Fragments Per Kilobase of transcript per Million fragments mapped)分析基因表达情况[10]。采用Ballgown对长瓣兜兰花蕾和花差异表达基因进行分析[10]。差异基因筛选条件设置为|log 2 foldchange|>1且校正P-value(FDR)<0.01。

2 结果与分析

2.1 原始测序数据质控

通过Illumina 测序平台测序后，分别获得484 702 56条和448 037 28条clean reads。两个样本的总的碱基数分别为665 508 257 3和619 826 070 2，GC含量占比分别为52.38%和52.64%， Q 20占比分别为99.11%和99.06%，Q 30占比分别为96.98%和96.85%。在碱基的组成分布情况中 (图1)，横轴表示reads的位置，纵轴表示四种碱基分布的百分比。理论上，在测序过程中，A、T和G、C的含量在每个循环上应分别相等，保持稳定不变，呈水平线。但是，由于通过随机引物扩增存在偏差的原因，会引起在测序得到的每个read前6～7个碱基有较大波动，但这种波动属正常情况。在碱基质量分布情况中 (图2)，横轴表示reads碱基所处位置，纵轴表示所有reads在该碱基位置上的质控得分 (0～40)。其中红色表示中位数，黄色代表25%～75%区间，触须代表10%～90%区间，蓝线处为平均数。高通量测序是双端测序，每条read的长度为150 bp。随着测序的进行会导致酶的活性逐步降低，因此，当测序达一定长度后，碱基的质量值会随之降低。但是，在本研究中，所有 reads在该碱基上的质控得分均处于绿色区域，表明数据质量可接受。

2.2 长瓣兜兰花器官转录组测序与拼接组装

由于长瓣兜兰基因组测序尚未完成，转录组测序数据只能进行无参考基因组的分析。下机数据获得后，进行组装测序数据，产生重组群和单一序列，重头组装可为数据处理和基因生物学功能的分析奠定基础。通过Trinity软件对clean data 进行从头组装，共计得到 170 807条transcript，由于一个基因可能会有两个或多个转录本，因此取最长的 transcript 作为unigene，所得的unigene总共为96 659 个。其中最长的转录本为 10 567/10 567 bp；最短的转录本201/201 bp；平均长度为 546.83/638.86 bp；组装转录本按从小到大的顺序排序，当转录本累加的长度占总长度的50%时，所对应的unigene和transcript长度分别为 923 bp、772 bp；当累加转录本的长度达到总长度的90%时， 所对应的unigene和transcript长度分别为241 bp 、273 bp (表1)。如图3所示，片段长度主要集中在 200～1 000 bp之间，其中转录本长度在 200～300 bp之间数量最多，其次是300～400 bp，长度越长所占比例越小。

表1 转录组组装数据长度分布

2.3 Unigene的功能注释

为了探究转录本所具有的生物学功能，对转录本序列进行基因功能注释和分类。使用 BLAST软件将 Unigenes 序列与NR、Swissprot、KEGG以及KOG数据库比对，获得Unigenes 的注释信息。长瓣兜兰花器官最终获得注释信息的 Unigenes 有 61 629条，占总Unigene总数的64.43%。长瓣兜兰花器官转录组Unigene 在 NR、KOG、Swissprot、KEGG 数据库中被注释的基因分别占总 Unigene的 47.64%、29.69%、54.51%、5.12%。长瓣兜兰转录组 Unigene 在 CDD、KOG、NR、NT、PFAM、Swissprot、TrEMBL、GO、KEGG等数据库中被注释的基因数目分别为 33 589、28 405、45 568、56 635、23 870、52 141、44 973、54 934、4 893。

从图4可以看出， 长瓣兜兰转录组测序组装的 Unigenes 与其他物种 的 Unigenes 相似数量最多，达到了29.42%；其次是油棕(Elaeisguineensis)，占 NR 注释总序列的 20.54%； 海枣(Phoenixdactylifera)为第三。

在GO 注释分析中，Unigene 序列主要分为分子功能、细胞组分、生物过程 3 个大类。从图5可以看出，在所有差异表达基因中，参与分子功能的基因有22类，细胞组分的为22个，还有27类聚集在生物学过程中。其中，结合、催化活性、转运活性、 细胞过程、代谢过程、细胞 及细胞部分关联的Unigene 较多。

为进一步评价转录组的完整性和注释的有效性，对 Unigene 进行KOG分类。共发现有 28 405 条 Unigen 得到注释，获得功能分类 25 个，在这些分类中，以信号转导机制类 (4 252 条)最多，一般功能预测类 (3 575条)细胞运动最少 (31个) (图6)。结果说明，在长瓣兜兰花发育的过程中通过大量的生物信号转导来完成其生物学过程。

为了进一步了解基因的生物学功能，采用KEGG数据库对长瓣兜兰可能参与的生理生化反应途径进行预测。分析结果表明，4 893 条 Unigene 得到注释，共参与了 23 类327个代谢途径 (图7)，其中包含Unigene 较多的途径有：碳水化合物代谢途径、能量代谢途径、脂质代谢途径等。代谢通路富集表明在长瓣兜兰花发育过程中，碳水化合物、能量和脂质代谢旺盛。

2.4 基因表达分析

通过FPKM分析，长瓣兜兰花蕾和花朵均表达的基因有 61 982个，只在其花蕾期特异表达的基因有27 488 个，只在花朵中特异表达的基因有6 000个。根据差异基因筛选条件，在 2 个发育时期筛选出 6 967 个差异表达基因，其中花朵相对花蕾上调表达的有2 417个，下调表达的有4 550个。在差异表达的基因中，存在大量与花发育相关的基因，但部分基因在花蕾和花朵的表达中存在显著差异，比如，相对于长瓣兜兰花蕾期，AG基因和C2H2-ZFP基因在花朵期表达量显著下调，而AGL80、WUS，FPA基因和FLC基因的表达量显著升高(图8)。

2.5 SSR标记开发

通过用MISA 软件对 Unigene 进行 SSR分析，发现在30 709 条Unigene中共有 7 613 条有 SSR，包括了 8 160个由一至六个核苷酸重复序列组成的 SSR位点。SSR 丰富度依次为二核苷酸 (2 567,33.72%)、一核苷酸 (2 445,32.12%)、三核苷酸 (1 988,26.11%)、四核苷酸 (72,0.95%)、六核苷酸 (19,0.25%)和五核苷酸 (16,0.21%)。长瓣兜兰SSR特征分析为开展长瓣兜兰及兜兰属植物分子标记及遗传图谱构建奠定基础。

3 讨 论

采用RNA-seq技术对长瓣兜兰两个不同发育时期的花器官进行了转录组无参基因组分析，为解析长瓣兜兰花发育及开花的分子机制奠定基础，也为利用分子生物学手段进行长瓣兜兰的花期调控和分子育种提供可利用的基因资源。

通过GO功能预测，发现共有54 934个Unigene得到注释，通过分类可将其分为分子功能、生物学过程和细胞组分三个类别，共71个小组，其中，结合、催化活性、转运活性、细胞过程、代谢过程、细胞及细胞部分关联的Unigene较多，这可能是由于在长瓣兜兰花芽分化进程中细胞不断增殖导致花芽内代谢活动旺盛引起的。这与橡胶树花序、雄蕊及雌蕊的注释结果大致相同[11]，说明在细胞代谢活动比较旺盛的生殖器官发育时期，大量的功能基因处于表达状态。

在差异表达的基因中，与花发育相关的基因主要有AG、AGL80、C2H2-ZFP、WUS、FPA以及FLC。在这些基因中，相对于长瓣兜兰花蕾期，AG和C2H2-ZFP在花朵期表达量显著下调，而AGL80、WUS，FPA和FLC的表达量显著升高。花朵时期，大量的基因表达下调，可能是由于到花朵时期，花的各个部分器官分化已完成，调控器官分化的基因已完成使命导致表达量下降。有研究表明，AG属于ABC模型中C功能基因，在调节花朵个体发育的后半部分起到至关重要的作用，表现在调控花柱和心皮的发育以及花朵发育的终止上[12]；而C2H2-ZFP可通过影响细胞增殖和分裂过程中的激素信号，直接或间接调控花器官发育中的ABC模型功能基因[13]；WUS则对花分生组织的维持和确定性至关重要[12]；FPA蛋白组主要控制开花时间；而FLC能抑制春化期花期的起始。在本研究中得到与花发育相关的基因在花发育过程中具体表达量如何，是如何起作用的，将在后期实验中进行研究，以更好地发掘出长瓣兜兰花发育和开花的相关基因。

优良品种的选育离不开分子标记的手段，开发SSR分子标记对品种的选育具有重要意义。通过表型数据的测定和观察作为育种的主要参考依据，用分子标记进行辅助选育，可以培育出品质更优、表型更好的优良品种[14]。在花卉植物中，已经有很多植物进行了SSR分子标记的开发，如朱顶红[15]、百合[16]、兰科植物寒兰[17]、亨利兜兰[18]以及同色兜兰(Paphiopedilumconcolor)均进行了分子标记的开发，这对育种具有极大的意义。在本研究中，采用软件 MISA(1.0 版)对 Unigene 进行 SSR 检测，发现7 613 条有 SSR，共搜索到8 160个SSR 位点，这比杂交兰[19]中的SSR位点少得多。本研究虽通过软件预测了SSR分子标记，但在位点多态性方面可能还会存在差异。因此，后期还需结合实验进一步验证SSR分子标记的准确性。

4 结 论

本研究用 Illumina Hiseq测序平台对长瓣兜兰花器官进行无参转录组分析，共获得95 659个 Unigene，其中，有45 568 个 Unigene 在 NR数据库中获得注释, 其中有9 345个Unigene注释到油棕中；在GO数据库中注释到54 934个Unigene，可将其分为分子功能、生物学过程和细胞组分三大类和71个亚类；在KOG 数据库中有28 405个 Unigene 获得注释，可分为 25 个功能区域；KEGG 数据库中有注释到的4 893 条 Unigene 参与了 23 类327个KEGG代谢途径，可以全面了解长瓣兜兰花器官的代谢途径信息；分析了在花蕾期和花朵期的差异表达基因，共有6 967 个差异表达基因，其中花朵相对花蕾上调表达的有2 417个，下调表达的有4 550个；在长瓣兜兰花器官转录组中存在8 160 个 SSR 位点。本研究通过转录组测序，获得了大量长瓣兜兰花器官基因序列信息，初步了解长瓣兜兰花器官发育过程中基因的表达情况，这为后续进一步深入开展长瓣兜兰基因功能研究、花发育的分子机制及 SSR 分子标记开发奠定了基础。