药用植物猴耳环的转录组测序分析*

谭志强 黄炬峰 李福明 邓乐平 毛积鹏,2

（1.台山市红岭种子园，广东 江门529223；2. 华南农业大学 林学与风景园林学院/广东省森林植物种质创新与利用重点实验室，广东 广州 510642)

猴耳环Pithecellobium clypearia是我国南方一种重要的药用植物，具有抗菌、抗炎、抗氧化、抗流感、抗肿瘤、降血糖、降血脂和免疫等作用[1-8]。特别是其高效的抗炎作用在临床上被广泛应用于呼吸道感染、咽喉炎、扁桃体炎和肠胃炎等炎症的治疗。目前已研发出了猴耳环消炎片和猴耳环消炎胶囊等产品[9-11]。猴耳环的主要药用部位为干燥的嫩枝和叶。近年来相关学者利用高效液相色谱（HPLC）和色谱—质谱联用（LC-MS）等技术对猴耳环的主要活性成分进行了测定。结果表明：猴耳环的主要活性成分为没食子酸、槲皮素、槲皮苷、杨梅苷、吡喃鼠李糖苷和表没食子儿茶素等黄酮类化合物[12]。虽然到目前为止已经开展了许多关于猴耳环的研究，但主要集中于其生长特性研究、繁殖技术优化和活性成分的鉴定等方面[13-15]，其转录组测序信息尚未见报道。

转录组测序能在缺乏基因组信息的条件下，有效地挖掘功能基因并揭示其生物学特性和基因的内在关系。此外，转录组测序获得的EST-SSR分子标记也是开展群体遗传多样性分析和分子育种的重要基础。目前，如凤丹Paeonia suffruticosacv. Feng Dan[16]、半枫荷Semiliquidambar cathayensis[17]和罗布麻Apocynum venetum[18]等药用植物，大豆Glycine max[19]、野生大豆Glycine soja[20]和蒺藜苜蓿Medicago truncatula[21]等豆科植物已经完成了转录组测序分析，并鉴定出了一批和萜类、黄酮类和生物碱类等主要活性成分生物合成与代谢相关的候选基因。本研究利用Illumina HiSeq 4 000 高通量测序平台对猴耳环主要药用部位嫩枝和叶的混合样本进行转录组测序分析。旨在建立猴耳环的转录组数据库，为猴耳环群体遗传多样性分析、功能基因的挖掘和主要活性成分生物合成与代谢调控分子机制提供依据。

1 材料与方法

1.1 实验材料

猴耳环材料来自于广东省台山市红岭种子园。于2019 年9 月份取样，分别采取了2 a 生且生长旺盛的猴耳环植株的嫩枝和叶片组织，剪碎后放入50 mL 无菌的离心管中，随即放入液氮中，-80 ℃保存至RNA 的提取。

1.2 猴耳环总RNA 提取

猴耳环嫩枝和叶片样本总RNA 的提取根据自我改良的CTAB 法[22]进行，具体方法步骤如下：① 将约150 mg 的猴耳环样本组织放入已经灭菌的研钵中，在液氮下快速充分研磨成粉末状，随后立即转移到1.5 mL 灭菌的离心管中；② 往离心管中加入1 mL 自配的CTAB 裂解液（100 mM Tris-HCl，25 mM EDTA，2 M NaCl，2% CTAB，2% PVP，0.5 g/L 亚精胺，pH=8） 和20 μL 的β-巯基乙醇，充分混匀后在65 ℃的金属浴中温浴20 min，期间混匀2～3 次；③温浴后，4 ℃条件下，12 000 rpm 离心8 min，取上清液转移至另一个新的1.5 mL 离心管中（尽量避免吸取到沉淀物，可放弃一小部分上清不取）；④ 加入等体积的氯仿/异戊醇（24 : 1），充分混匀后，4 ℃条件下，12 000 rpm 离心5 min，收集上清液于另一个新的1.5 mL 离心管；⑤ 再次加入等体积的氯仿/异戊醇（24 : 1），充分混匀后，4 ℃条件下，12 000 rpm 离心5 min，收集上清液并加入1/4 体积的10 M LiCl 溶液，于4 ℃条件下静置5～6 h；⑥ 随后在4 ℃条件下，12 000 rpm 离心10 min，弃去上清液；⑦ 沉淀物先后用500 mL已提前-20 ℃预冷的75%乙醇和无水乙醇洗涤两次，4 ℃条件下，12 000 rpm 离心5 min，弃上清夜后置于室温条件下干燥5 min；⑧干燥后，加入30～50 μL 的DEPC 水溶解，随后进行总 RNA 浓度、纯度和完整性的检测，符合要求的总RNA 样品在-80 ℃条件下保存至cDNA 文库的构建。

1.3 cDNA 文库构建、测序及组装

cDNA 文库的建库流程和具体的操作步骤参考Foucart 等[23]的方法执行，使用广州瑞科基因科技有限公司（Science Corporation of Gene）的 Illumina HiSeq 4 000 测序平台对所构建的 cDNA 文库进行双末端测序。利用fastx_toolkit 软件中的fastx_clipper 工具和 fastq_quality_filter 工具去除测序接头盒低质量的reads 后，根据Hass 等[24]的 Trinity方法进行高质量 reads 的组装。

1.4 Unigenes 序列比对和功能注释

利用 BLASTn 软件对所组装获得的unigenes在Uniprot、Nr、KEGG、KOG 和GO 公共数据库中进行同源序列比对分析，其中 E 值设置为小于10-5。获得与序列相似度最高的蛋白质进行猴耳环 unigene 的功能注释。利用 Blast2GO 软件[25]将在Nr 数据库中获得注释的 Unigene 进行GO 注释,并用WEGO 软件[26]对其进行GO 分类功能统计。使用BLASTx 软件并结合KEGG 数据库进行unigenes 的KEGG 通路分析。

1.5 SSR 分析

简单重复序列（Simple Sequence Repeat, SSR）是一类由几个核苷酸为基本单元多次串联重复而形成测DNA 片段，分布于整个基因组中，数量大多态性丰富，具有广泛的应用性。本研究利用MISA 软件对猴耳环Unigene 的SSR 位点进行了检测。SSR 筛选标准是：单核苷酸重复10 次及以上，二核苷酸重复6 次及以上，三核苷酸重复5次及以上，四核苷酸重复5 次及以上，五核苷酸重复5 次及以上，六核苷酸重复5 次及以上。

2 结果与分析

2.1 转录组测序和组装分析

经质量评估和低质量测序数据筛除后，共获得24 748 936 个高质量的reads，碱基数为7 424 680 800，Q30 含量为93. 89%，GC 含量为40.17%。利用Trinity 软件对高质量测序数据进行拼接和组装，一共获得134 692 个转录本，转录本总长度为244 535 825 bp，最大为17 499 bp，最小为184 bp，平均长度为1 593 bp，N50 为2 509 bp。筛除相同转录本后共获得63 299 个Unigenes，平均序列长度为1 117 bp，序列长度分布如图1所示。

图1 猴耳环63 299 个Unigene 序列长度的分布Fig.1 Distribution of unigene lengths for Pithecellobium clypearia

2.2 Unigenes 功能注释分析

将获得的63 299 个Unigene 与Nr、UniProt、KEGG、KOG 和GO 共5 个公共数据库进行比对分析（E-value ＜ 10-5）。结果表明，26 101 个Unigene（41.23%）获得功能注释，其中25 591个Unigene 比对到Nr 数据库；25 983 个Unigene比对到UniProt 数据库；8 530 个Unigene 比对到KEGG 数据库；13 974 个Unigene 比对到KOG 数据库；19 299 个Unigene 比对到GO 数据库。根据Nr 数据库比对到的最相似的基因进行物种分布统计。结果发现3 383、3 211 和2 716 个Unigene 分别比对到木豆Cajanus cajan、大豆和狭叶羽扇豆Lupinus angustifolius的同源序列中（图2）。

图2 猴耳环转录组Unigenes 的Nr 数据库注释物种分布Fig.2 Nr database annotated species distribution of unigenes of transcriptome for Pithecellobium clypearia

2.3 Unigenes GO 功能分析

GO 数据库从生物学过程、细胞组分和分子功能3 个方面预测基因产物的功能。利用Blast2GO软件对在GO 数据库中注释到的Unigenes 进行功能分类分析。结果发现，19 299 个Unigene 被归类到36 个GO 功能亚类中，其中注释到生物过程大类中的Unigenes 数量最多，其次是分子功能和细胞组分大类（图3）。在生物过程大类中注释到的Unigenes 主要分布在代谢过程、细胞过程、单一生物过程和生物调节等功能亚类中。在分子功能大类中，绝大部分被注释到的Unigenes 分布在绑定和催化活性两个GO 功能亚类中。而在细胞组分大类中，细胞、细胞膜和细胞器部分3 个功能亚类中分布的Unigenes 显著高于其它GO 亚类。

2.4 KOG 功能分类分析

为进一步验证Unigenes 注释结果功能分类的可靠性，将所注释到的Unigenes 进行KOG 功能分类分析，各KOG 分类簇中Unigenes 的分布情况如图4 所示。结果表明：共13 975 个Unigene 被注释到25 个KOG 功能类别簇中。其中一般功能预测类别中分布的Unigenes 最多（3 759,26.8%），其次是翻译后修饰、蛋白转运、分子信号和信号转导机制等KOG 功能类别。而在细胞核结构、细胞外结构和细胞运动3 个KOG 功能类别中分布的Unigenes 分别只有6、1 和1 个。

图3 猴耳环转录组Unigene 的GO 功能分类统计Fig.3 GO functional classification of unigenes of transcriptome for Pithecellobium clypearia

图4 猴耳环转录组Unigenes 的KOG 功能分布Fig.4 KOG functional annotation distribution of unigenes of transcriptome for Pithecellobium clypearia

2.5 KEGG 代谢通路分析

在KEGG 数据库中注释到的8 530 个Unigenes主要分布在231 个KEGG 代谢通路中。主要涉及到了代谢、遗传信息处理、环境信息处理、细胞过程及生物系统5 个大类通路和26 个子类通路中。结果表明，在KEGG 数据库中被注释到的Unigenes 主要分布在代谢通路大类中，而生物系统通路大类中分布的Unigene 最少。萜类化合物和黄酮类化合物为猴耳环的主要活性成分，分别通过萜类骨架和苯丙烷途径合成。本研究基于萜类和黄酮类化合物生物合成KEGG 通路及功能注释结果，鉴定到192 个和黄酮类化合物生物合成相关的Unigenes，主要定位在苯丙烷生物合成通路、黄酮类化合物生物合成通路和异黄酮生物合成通路中（表1）。共鉴定到136 个和萜类化合物生物合成相关的Unigenes，其中参与萜类骨架生物合成的Unigene 53 个、参与二萜类化合物生物合成通路的Unigenes 23 个，参与四萜类化合物类胡萝卜素生物合成通路的Unigene 29 个、参与倍半萜和单萜类化合物生物合成通路的Unigene 28 个（表1）。

2.6 SSR 位点分析

本研究利用MISA 软件对组装得到的134 692个转录本中的SSR 位点进行检测。结果在15 740个Unigenes 中共检测到45 573 个SSR 位点，发生频率为33.84%。其中5 928 个Unigenes 含有2 个及2 个以上的SSR 位点。在猴耳环转录本SSR 类型中共检测到单核苷酸至六核苷酸6 种重复类型。主要以单核苷酸重复（26 377, 57.87%）、二核苷酸重复（11 755，25.79%）和三核苷酸重复（6 853，15.04%）为主。四核苷酸重复、五核苷酸重复和六核苷酸重复只占总SSR 位点的1.29%。在单核苷酸重复类别中，出现频率最高的为A/T，占96.47%。出现频率较高的AT/TA、AG/TC 和CT/GA 三种重复类型分别占二核苷酸SSR 位点总数的26.47%、26.27%和23.23%。三核苷酸的主要重复类型为AAG/TTC、CTT/GAA 和AGA/TCT。四核苷酸中，TTTA/AAAT 和AAAG/TTTC 为其优势重复类型。然而在五核苷酸和六核苷酸SSR 位点中没有具显著优势的重复类型。

3 讨论与结论

随着测序技术的不断提高和测序成本下降，高通量转录组测序在无参考基因组物种功能基因挖掘和分子标记开发等方面的应用越发广泛。本研究利用RNA-seq 技术首次对猴耳环的主要药用部位嫩枝和叶片的混合样本进行了转录组测序分析。共获得2 474 893 个高质量reads，Q30 高达93.89%，表明猴耳环的cDNA 文库构建质量高。拼接和组装后共获得63 299 个Unigenes，平均序列长度为1 117 bp。与其它已完成转录组测序的药用植物相比，组装到的Unigenes 与凤丹[16]、五指毛桃Ficus hirta[27]、半枫荷[17]和罗布麻[18]等药用植物相近，初步说明猴耳环的转录组序列的拼接组装效果较好，测序质量可靠。但明显低于同为豆科的大豆[19]、野生大豆[20]和蒺藜苜蓿[21]等植物。这可能是因为猴耳环缺少基因组信息导致的。

表1 猴耳环黄酮类和萜类化合物生物合成途径相关基因Table 1 Related genes of flavonoids and terpenoids biosynthesis in Pithecellobium clypearia

表 2 猴耳环Unigenes 中SSR 重复单元的分布特征统计Table 2 Statistical information of the SSR motifs distribution in Pithecellobium clypearia

将猴耳环组装到的Unigenes 序列与5 个公共数据库进行比对分析。结果发现只有26 101 个Unigenes（41.23%）有功能注释，58.77%的Unigenes 没有获得功能注解信息。这与野三七Panax stipuleanatus[28]、 半夏Pinellia ternata[29]和罗布麻[21]等药用植物的转录组数据的功能注释结果相似。可能是由于短序列占比大、保守的核心序列信息不完整等因素导致的。成功注释到Nr 数据库中的25 591 个Unigenes，65%比对上了木豆、大豆和蒺藜苜蓿的等11 种豆科植物的同源序列。这可能是由于猴耳环同为豆科植物与其具有较近的亲缘关系，亦可能是因为木豆、大豆和蒺藜苜蓿等物种具有参考基因组信息。在GO、KOG 和KEGG 数据库中也分别注释到19 299、13 974 和8 530 个Unigenes，主要分布在一般功能预测和代谢过程及通路中，说明猴耳环主要药用部位嫩枝和叶组织中参与次生代谢产物生物合成与调控的基因丰富。萜类和黄酮类化合物为猴耳环的主要活性物质，本研究中共鉴定到192 和黄酮类化合物生物合成相关的Unigenes，主要参与苯丙烷、黄酮、黄酮醇、异黄酮和花青素等生物合成通路。136 个和萜类化合生物合成相关的Unigenes，主要参与萜类骨架、单萜、倍半萜、二萜和四萜化合物合成通路。这为揭示猴耳环萜类和黄酮类化合物合成及代谢调控网络提供了理论基础。SSR 标记因其通量大、操作简单和重复性好等特点已被广泛应用于铁皮石斛Dendrobium officinale[30]、木麻黄Casuarina equisetifolia[31]、枸杞Lycium barbarum[32]和黄芩Scutellaria baicalensis[33]等药用植物分子标记辅助育种和遗传图谱的构建。本研究利用MISA 软件在15 740 个Unigenes 中共检测到45 573 个SSR 位点，发生频率高达33.84%，显著高于连翘Forsythia suspensa[34]、人参Panax ginseng[35]、鱼腥草Houttuynia cordata[36]和杜仲Eucommia ulmoides[37]等药用植物中SSR 出现的频率。猴耳环转录组中不同重复类型的SSR 数量具有较大差异，主要以单核苷酸、二核苷酸和三核苷酸为主。在单核苷酸中，大部分为A/T 型SSR 位点，G/C 含量少。导致这种差异的原因可能是因为甲基化C 残基变为T 碱基。此外二核苷酸和三核苷酸SSR 重复类型也占据了很大比例，这与罗布麻和杜仲等药用植物转录组数据的SSR 重复类型比例相类似。本研究鉴定到SSR 分子标记可为分析猴耳环的遗传多样性、构建遗传图谱、挖掘功能基因和分子标记辅助育种等研究提供理论基础。