杜梨叶片转录组分析

陈 燕， 郭 聪， 王 莹， 李玉娟

(江苏沿江地区农业科学研究所，江苏南通 226541)

杜梨(Bunge)隶属于蔷薇科梨属，原产于我国，耐寒耐旱，是常用的梨砧木。杜梨的枝叶、树皮、根和果实均有药效，能消食止痢，润肺止咳。杜梨作为我国本土树种，耐盐碱性强，具有林用、观赏、果用等许多优良特性，是一个值得深度开发的树种。将杜梨种植在难以治理的盐碱化等地区进行绿化，可以缓解环境污染、水土流失，扩大绿化面积，提高观赏性和生态效益。近年来，江苏省极力推进沿海大开发，沿海滩涂绿化与生态保护建设方兴未艾。

国内外对杜梨的生物学研究多集中于耐盐、耐旱的分子机制上，与叶色相关的研究只有马春晖等报道过，即以我国梨属主要植物种(杜梨、豆梨、西洋梨)以及近缘属植物榅桲为试材，通过植物叶片色泽和类别进行相关性分析，发现它们之间存在一定的相关性。然而，尚未见关于其叶片转录组测序方面的研究报道。

RAN-Seq测序技术能够分析植物整体转录活动，是研究彩叶植物呈色机制的重要方法。郭聪等发现了7个转录因子通过调控花青素生物合成途径来改变美国红枫叶色。Gang等通过转录组测序技术发现黄绿叶突变体桦树中光合作用天线蛋白代表最重要的富集途径，其表达水平对叶片颜色的形成至关重要。褚江涛对三角枫新品种齐鲁金叶片进行了转录组分析，共鉴定出75个差异表达基因，并通过转录因子预测，筛选出180个下调的转录因子基因。陆小雨等对红花槭特殊单株不同色叶进行代谢组和转录组分析，发现秋季花青素苷通路中大量基因的表达量上调是红花槭叶片变色的主要驱动因子。申建双研究了金叶连翘叶片的呈色机制，对转录组数据分析发现，与金叶连翘叶色形成及响应光照度有关的基因45个、转录因子23个。其中，15个基因表达不受光照度影响，30个基因和23个转录因子的表达均受光照度的影响。最终筛选出6个基因和1个转录因子在黄叶型连翘中显著低表达或不表达，为连翘叶色变黄及响应不同光照度的关键候选基因。

综上，转录组测序已广泛应用于植物叶片呈色机制的分子研究中，但目前还未见杜梨叶色转录组测序分析相关的报道。本研究以转色期杜梨青叶、花叶、红叶为试验材料，分析比较基因表达差异，以期开发出可用于耐盐地开发，集绿化、林用等多功能的红叶杜梨园林植物新品种(系)，做到春天观花、夏天观果、秋后赏红叶，为城市美化与江苏沿海地区大开发生态环境建设作出贡献。

1 材料与方法

1.1 试验材料

2019年11月下旬于江苏沿江地区农业科学研究所杜梨培育基地采集杜梨单株上不同分枝的红叶、花叶、青叶为试验材料(图1)，每组叶片分别取3次重复，编号见表1；将叶片用锡箔纸包裹好后，迅速置于液氮中。

杜梨高通量测序和分析由北京百迈克生物科技有限公司负责进行。

1.2 试验方法

1.2.1 杜梨总RAN提取及质量检测 采用试剂盒法提取RNA并采用超微量分光光度计对样品RNA质量、浓度、完整性进行检测。杜梨叶片的总 RNA 的提取是采用试剂盒NanoDrop 2000(Thermo Fisher Scientific，美国)。RNA 样品的完整性是通过 NanoDrop6000 (AgilentTechnologies，美国)进行检测。

1.2.2 cDNA文库构建 以合格的mRNA为模板，通过六碱基随机引物、AMPure XP beads核酸纯化试剂盒和 PCR 富集等相关工具和技术获得 cDNA 文库。

1.2.3 文库质控及测序 使用Q-PCR方法对cDNA文库的有效浓度(文库有效浓度>2nmol/L)进行准确定量；库检合格后，用Illumina平台进行测序。

2 结果与分析

2.1 杜梨转录组测序与组装

对9个样品的转录组测序和分析，共获得67.78 Gb测序数据，各样品总碱基数均达到6.79 Gb，测序数据质量大于或等于30的碱基所占的百分比(Q30)在91.73%及以上，GC含量百分比均在46.22%及以上(表1)，可用于下一步组装。利用StringTie软件对比对到参考基因组上的片段(mapped reads)进行拼接组装，共得到60 611个基因，并与原有的基因组注释信息进行比较，经筛选共发掘出1 059个新基因。

表1 杜梨样品测序数据评估

2.2 差异表达分析

皮尔逊相关系数(Pearson’s correlation coefficient，)是评估生物学重复相关性的重要指标，越接近1，说明相关性越强。从图2可以看出，CK1、CK2、CK3，B1、B2、B3，H1、H2、H3之间的均大于0.900，表明样本间相关性强，所有样本数据均可用于后续分析。

将(fold change，差异倍数)≥2且(false discovery rate，错误发现率)<0.01作为筛选标准，对杜梨3组不同颜色叶片的差异表达基因进行对比分析。从表2、图3可以看出， 青叶和花叶共有57个DEGs表达量发生变化，其中15个DEGs上调，42个DEGs下调；青叶和红叶共有2 084个DEGs表达量发生变化，其中1 447个DEGs上调，637个DEGs下调；红叶和花叶对比下，共有1 795个DEGs表达量发生变化，其中568个DEGs上调，1 227个DEGs下调。

表2 杜梨差异表达基因数量统计

将每组差异基因进行比对，通过韦恩图来了解差异基因在每组中的分布情况(图4)，青叶与花叶、青叶与红叶之间共有20个差异表达基因；红叶与青叶、花叶之间共有1 258个差异表达基因，花叶与青叶、红叶之间共有20个差异表达基因，三者共有9个差异表达基因。

2.3 基因功能注释

将全部基因序列与COG(clusters of orthologous groups，直系同源群)、GO(gene ontology，基因本体)、KEGG(kyoto encyclopedia of genes and genomes，京都基因与基因组百科全书)、KOG(eukaryotic orthologous groups，真核生物同源组)、Pfam(protein family，蛋白家族)、Swiss-Prot(manually annotated and reviewed protein sequence database)、eggNOG(evolutionary genealogy of genes：non-supervised orthologous groups)、NR(non-redundant，非冗余数据库)数据库进行比对(表3)发现，共57 669个基因得到功能注释，注释率为95.15%。其中，序列长度在≥300～1 000 bp之间的基因占总注释基因的37.91%(21 862/57 669)，序列长度≥1 000 bp的占59.71%(34 436/57 669)。对差异表达基因功能注释进行统计，发现青叶与花叶、青叶与红叶、花叶与红叶分别有56、2 064、1 769个差异表达基因被注释，分别占98.25%(56/57)、99.04%(2 064/2 084)、98.55%(1 769/1 795)，均超过98%(表4)，表明差异表达基因的高注释率可能为杜梨色叶研究和新品种选育提供了理论基础。

表3 unigene功能注释统计

表4 差异表达基因功能注释统计

2.4 差异表达基因GO分类和功能富集

将杜梨全部基因序列与GO 数据库比对，发现30 403个基因得到了功能注释，其中差异表达基因得到功能注释数最多的为青叶与红叶，共1 394个；注释率最高的为青叶与花叶，为76.79%(表5)。注释的基因涉及49个功能组，主要集中在膜(细胞组分)、催化活性(分子功能)、代谢过程(生物学过程)等(图5)。通过GO分类，把杜梨基因按照功能进行分类，为杜梨抗逆性等研究提供了依据。

表5 GO数据库三大功能主类DEG数

2.5 差异表达基因KEGG注释和富集分析

通过KEGG数据库比对，共有16 594个杜梨基因得到注释。其中，青叶和花叶有11个DEGs可以分配到16个KEGG通路中，6个差异表达基因呈上调趋势，5个基因呈下调趋势，显著富集在碳代谢、内质网中的蛋白质加工、植物-病原体相互作用以及倍半萜和三萜生物合成中；青叶和红叶共有383个差异表达基因得到注释信息，分配到105个KEGG通路中，其中262个差异表达基因呈上调趋势，121呈下调趋势，显著富集在植物-病原体相互作用、碳代谢等通路中；红叶和花叶共有318个差异表达基因得到注释信息，分配到106个KEGG通路中，其中119个差异表达基因呈上调趋势，199个呈下调趋势，显著富集在植物-病原体相互作用、碳代谢等通路中(图6、图7)。综上，碳代谢、植物-病原体相互作用一直是杜梨叶片转色过程中排名前三的代谢通路，可能参与杜梨叶片转色。

2.6 与叶色相关基因分析

类黄酮是植物中通过苯丙烷途径合成的次生代谢产物，其产生3种主要化合物：黄酮醇，原花青素和花青素，与叶色息息相关。在杜梨叶片转色期间，从红叶、花叶、青叶等3种不同叶色的叶片中筛选出10个与类黄酮相关的基因，分别属于UDPGT(Chr2.g42475、Chr11.g09871、Chr11.g10237、Chr14.g49947、Chr15.g03978)、NAD dependent epimerase/dehydratase family(Chr5.g05869)、2OG-Fe(Ⅱ) oxygenase superfamily(Chr5.g09527)、ABC transporters(Chr12.g36590)、Transferase(Chr17.g27158)、Chalcone(Scaffold29.g59180)等6个基因家族，主要通过类黄酮的生物合成过程、花青素还原酶活性调节、黄酮醇合酶活性调节、花青素在紫外线下于组织中积累、黄酮类生物合成过程等机制对杜梨叶色进行调控(表6)。

表6 与叶色相关差异表达候选基因及GO注释

为了探索这10个基因在杜梨叶片转色期间的表达模式，本研究对这10个基因的表达水平进行分析，结果(图8)发现，每个基因都呈现出不同的表达水平，大部分转录因子在红叶中的表达量大于在青叶和花叶中的表达量。

2.7 SNP分析

本研究对所筛选出的SNP位点数目、转换类型比例、颠换类型比例进行统计分析，从表7可以看出，共筛选出452 501个SNP位点，发生频率为1/411。其中，碱基转换位点共277 789个，占SNP总数61.39%；碱基颠换位点共174 712个，占38.61%。所有的变异类型中，转换类型远远高于颠换类型，本研究结果极大地丰富了杜梨SNP位点信息，为今后杜梨种质资源筛选、新品种选育等研究提供了理论基础。

表7 SNP类型分析

3 讨论与结论

近年来，RNA-Seq、Chip-Seq等新一代技术取得了前所未有的发展，已经通过转录组、蛋白组、代谢组等不同层面和角度对植物呈色机制进行深度解析，为彩叶植物新品种培育等研究奠定了基础。其中，紫薇、美国红枫、银杏、花叶矢竹、墨兰、凤香果等彩叶植物已进行了转录组测序。本研究对杜梨叶片(青叶、花叶、红叶)进行了转录组分析，共获得了67.78GB总碱基，Q30碱基比例均大于91.73%。通过组装得到60 611个基因，包括1 059个新基因，为今后杜梨叶片呈色机制的研究、相关基因的挖掘以及新品种的选育提供了理论基础。

相关研究表明，同一个转录分析中unigene的序列越长，匹配率越高。本研究共发现57 669个基因得到功能注释(注释效率为95.15%)，其中有21 862个已注释基因的序列长度在≥300～1 000 bp 之间，34 436个已注释基因序列长度≥1 000 bp，与前者研究相似。类黄酮物质在植物叶片变色期起到重要的作用，其包含的一类物质——花色苷是植物叶片呈色的主要因素，其中UDPGT家族在花色苷合成途径中起重要的作用。本研究共筛选出10个与类黄酮相关的转录因子，其中属于UDPGT家族的共有5个转录因子，推测这5个UDPGT转录因子可能与郭聪等在美国红枫转录组中和王莹等在紫薇转录组中筛选的UDPGT转录因子具有同源性，在杜梨叶片变色期起到重要的作用。另外，在杜梨转录组中发现Chr5.g09527在3种叶片中的表达量显著高于其他的转录因子，但在杜梨叶片转色期的作用还需进一步分析。

SNP分子标记在品种鉴定、功能基因定位等方面具有重要的作用。本研究从杜梨叶片转录组中共搜索到452 501个SNP位点，发生频率为1/411，低于水稻(1/89)、玉米(1/61)等模式植物以及梨(1/344)、苹果(1/288)、葡萄(1/64)、火龙果(1/204)等果树，可能与基因组测序深度深、物种间的遗传背景差异大有关，也可能与SNP检测方法有关。另外，在杜梨叶片转录组中的所有SNP变异类型中，碱基转换类型(61.39%)高于颠换类型(38.61%)，这与梨(31.50%、31.66%)等大多数植物研究结果相似。碱基转换类型比例高可能与甲基化等原因相关。

本研究对杜梨转色期叶片进行了转录组测序和分析，共获得了67.78 Gb测序数据，通过组装共获得60 611个基因，其中包含1 059个新基因；使用Blast软件将全部基因序列与8个数据库进行比对，统计到57 669个基因得到功能注释；GO分析显示30 403个基因得到功能注释，共涉及49个功能组，主要集中在膜(细胞组分)、催化活性(分子功能)、代谢过程(生物学过程)等中；KEGG富集分析发现16 594个基因得到注释，要对碳代谢和植物-病原体相互作用等生物学过程进行调节；在杜梨不同颜色的叶片中共筛选出10个与类黄酮相关的转录因子，主要通过类黄酮的生物合成过程、花青素还原酶活性调节、黄酮醇合酶活性调节、花青素在紫外线下于组织中积累、黄酮类生物合成过程等机制对杜梨叶色进行调控，并且大部分转录因子在红叶中的表达量大于在青叶和花叶中的表达量；在杜梨叶片转录组中共搜索到452 501个SNP位点，发生频率为1/411，所有的变异类型中转换类型远远高于颠换类型。