药用植物粗茎鳞毛蕨转录组测序分析
黄育通1, 邢艳萍1, 门文效1, 黄彦昌1, 杨燕云1, 康廷国1, 窦德强1,2, 许亮1,2     
1. 辽宁中医药大学药学院, 辽宁大连 116600;
2. 辽宁省中药质量及资源开发专业技术创新中心, 辽宁沈阳 110847
摘要: 为获取粗茎鳞毛蕨(Dryopteris crassirhizoma)的转录组信息并分析其基因表达特征,本研究采用Illumina HiSeqTM 2000 150PE高通量测序平台, 对其叶、叶柄、叶轴、根和根状茎进行转录组测序分析。结果表明:经Trinity拼接并去冗余后,共获得147 690个unigenes,在核酸序列库(Nucleotide sequence database,NT)、非冗余蛋白质数据库(Non-Redundant protein database,NR)、Swiss-Prot蛋白质序列数据库(Swiss-Prot protein sequence database,Swiss-Prot)、京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes,KEGG)、直系同源蛋白质群数据库(Clusters of Orthologous Groups of proteins,COG)和基因本体库(Gene Ontology,GO)中获得注释的unigenes数量分别为91 758、57 578、60 535、56 670、41 446、63 027个。GO分类注释得到3大类56个分支,KEGG代谢通路注释到129条标准通路,其中23条为次生代谢通路。蛋白质编码框序列94 731个,以短序列为主。简单重复序列(Simple Sequence Repeat,SSR)分析发现28 079个SSRs,以二核苷酸SSRs最多,占总SSRs数量的70.12%。单核苷酸多态性(Single Nucleotide Polymorphism,SNP)分析检测到叶、叶柄、叶轴、根和根状茎的SNP数量分别为23 139、21 246、28 280、55 409、17 084个。叶柄与根状茎差异基因表达量分析显示上调基因13 829个、下调基因24 250个。GO功能富集分析结果表明,在叶柄和根状茎组织中共富集到226个与苯丙烷类代谢途径相关的差异表达基因。本研究结果可为后续深入探究粗茎鳞毛蕨次生代谢通路、间苯三酚类与黄酮类等物质的生物合成及基因功能提供转录组数据支持。
关键词: 粗茎鳞毛蕨    转录组    功能注释    代谢通路    差异基因    
Transcriptome Sequencing Analysis of the Medicinal Plant Dryopteris crassirhizoma
HUANG Yutong1, XING Yanping1, MEN Wenxiao1, HUANG Yanchang1, YANG Yanyun1, KANG Tingguo1, DOU Deqiang1,2, XU Liang1,2     
1. School of Pharmacy, Liaoning University of Traditional Chinese Medicine, Dalian, Liaoning, 116600, China;
2. Liaoning Province Technical Innovation Center for Traditional Chinese Medicine Quality and Resource Development, Shenyang, Liaoning, 110847, China
Abstract: To obtain transcriptomic information of Dryopteris crassirhizoma and analyze its gene expression characteristics, this study performed transcriptome sequencing on leaves, petioles, leaf axes, roots, and rhizomes by using the Illumina HiSeqTM 2000 150PE high-throughput sequencing platform.The results showed that after assembly with Trinity and removal of redundancy, a total of 147 690 unigenes were obtained.Annotation was performed against the Nucleotide sequence database (NT), Non-Redundant protein database (NR), Swiss-Prot protein sequence database (Swiss-Prot), Kyoto Encyclopedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups of proteins (COG), and Gene Ontology (GO), with 91 758, 57 578, 60 535, 56 670, 41 446, and 63 027 unigenes annotated, respectively.GO annotation categorized the genes into three major categories comprising 56 subcategories.KEGG pathway annotation identified 129 standard metabolic pathways, among which 23 were classified as secondary metabolic pathways.A total of 94 731 Coding Sequences (CDSs) were predicted, with the majority being short-length sequences.Simple Sequence Repeat (SSR) analysis identified 28 079 SSRs, with dinucleotide SSRs being the most abundant, accounting for 70.12% of the total SSRs.Single Nucleotide Polymorphism (SNP) analysis detected 23 139, 21 246, 28 280, 55 409, and 17 084 SNPs in the leaves, petioles, leaf axes, roots, and rhizomes, respectively.Differential gene expression analysis between petioles and rhizomes identified 13 829 upregulated and 24 250 downregulated genes.GO functional enrichment analysis revealed that 226 differentially expressed genes related to the phenylpropanoid biosynthesis pathway were enriched in both petioles and rhizomes.This study provides transcriptomic data support for future in-depth investigations into the secondary metabolic pathways of D.crassirhizoma, particularly the biosynthesis of phloroglucinols and flavonoids, and related gene functions.
Key words: Dryopteris crassirhizoma    transcriptome    functional annotation    metabolic pathways    differentially expressed genes    

粗茎鳞毛蕨(Dryopteris crassirhizoma)别名绵马鳞毛蕨、野鸡膀子,是一种药用蕨类植物,隶属于鳞毛蕨科(Dryopteridaceae)鳞毛蕨属(Dryopteris),主要分布在我国河北东部及东北地区[1]。贯众首载于《神农本草经》[2],是一种多基原中药材,其基原植物种类高达50余种[3],其中粗茎鳞毛蕨的干燥根状茎和叶柄残基被《中华人民共和国药典(一部)》(2020年版)[4]以绵马贯众为名收录。绵马贯众味苦性微寒,具有清热解毒、驱虫等功效[4]。现代医学研究表明,粗茎鳞毛蕨所含的间苯三酚类、黄酮类等生物活性物质具有抗病毒、抗肿瘤和抗炎等作用[5-7],其中间苯三酚类化合物为其特征性化合物[8]。绵马贯众在数次疫情防控期间均消耗量大,甚至一度脱销,是著名中成药连花清瘟胶囊的主要成分之一[9-10]。因此,绵马贯众的资源开发及其相关基础研究具有广阔的前景和重要的应用价值。

转录组的概念是由Charlles Auffary于1996年提出[11],与基因组相比,转录组具有时间性和空间性。转录组测序可揭示生物生长的内在规律,如次生代谢、转录调控等规律[12]。随着以寡核苷酸连接为代表的第二代测序技术和以纳米孔单分子为代表的第三代测序技术在近年来的迅猛发展[13],越来越多的药用植物,如杜仲(Eucommia ulmoides)[14]、朱砂根(Ardisia crenata)[15]、黄三七(Actaea vaginata)[16]、岩黄连(Corydalis saxicola)等[17],相继完成了转录组测序。然而,迄今尚未见有关粗茎鳞毛蕨转录组分析的报道,这不利于这一东北特色蕨类药用植物的进一步开发利用。因此,本研究采用Illumina HiSeqTM 2000 150PE高通量测序平台对粗茎鳞毛蕨叶、叶柄、叶轴、根和根状茎5种组织进行转录组测序,旨在揭示其转录组的整体表达特征,为后续基因挖掘及次生代谢途径解析提供数据支持。

1 材料与方法 1.1 材料

供试植物材料收集于2022年9月辽宁省鞍山市铁东区倪家台村(123°8′53.98″E,41°3′11.39″N),样品经辽宁中医药大学许亮教授鉴定为粗茎鳞毛蕨(Dryopteris crassirhizoma)。取同一样品植株的叶、叶柄、叶轴、根和根状茎5个部位新鲜组织,经液氮速冻后置-80 ℃冰箱保存备用。

1.2 RNA提取与文库构建

使用德国QIAGEN公司的RNeasy Mini Kit试剂盒分别提取粗茎鳞毛蕨叶、叶柄、叶轴、根和根状茎5种组织的总RNA。总RNA纯度使用纳米分光光度计(NanoPhotometer®,N60,德国Implen公司)检测,总RNA浓度使用美国赛默飞世尔科技公司的Qubit® RNA测定试剂盒和Qubit® 2.0荧光计(Life Technologies,美国CA公司)检测。总RNA完整性使用RNA 6000 Nano试剂盒(美国Agilent公司)评估。建库使用mRNA-seq v2建库试剂盒(南京诺唯赞生物科技股份有限公司),按试剂盒使用说明书进行,索引编码样品的聚类使用RNA Adapters set1/set2 for Illumina(南京诺唯赞生物科技股份有限公司)。

1.3 转录组测序与组装

于2022年12月采用Illumina HiSeqTM 2000 150PE高通量测序平台对粗茎鳞毛蕨5种组织进行转录组测序。每种组织设置3个生物学重复,每个重复样本均独立提取总RNA并建库进行测序,测序质量利用FastQC和Trimmomatic软件进行数据评估。测序原始数据上传至NCBI(National center for biotechnology information)数据库,并获得BioProject编号:PRJNA1302105(https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1302105)。测序结果使用短reads组装软件Trinity进行组装[18],组装得到的unigenes,首先使用Tgicl将其去除冗余和进一步拼接,再对这些序列进行同源转录本聚类,得到最终的unigenes。最后,将unigenes序列与核酸序列库(Nucleotide sequence database,NT)、非冗余蛋白质数据库(Non-redundant protein database,NR)、Swiss-Prot蛋白质序列数据库(Swiss-Prot protein sequence database,Swiss-Prot)、京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes,KEGG)和直系同源蛋白质群数据库(Clusters of orthologous groups of proteins,COG)做Blast比对分析(E值< 1×10-5),取比对结果最好的蛋白确定unigenes的序列方向。如果不同库之间的比对结果有矛盾,则按NT、NR、Swiss-Prot、KEGG和COG的优先级确定unigenes的序列方向,跟以上数据库皆比对不上的unigenes用ESTScan软件确定序列的方向[19-20]。能确定序列方向的unigenes给出其5′→3′方向的序列,无法确定序列方向的unigenes给出组装软件得到的序列。经过拼接、组装和聚类等分析后,得到unigenes和contig结果。

1.4 转录组功能注释

原始测序数据的生物信息分析参照从NCBI等数据库下载的参考序列[21]。通过Blast将unigenes序列与NT、NR、Swiss-Prot、KEGG和COG数据库进行比对分析(E值<1×10-5)。利用Blast2GO软件得到unigenes的基因本体库(Gene ontology,GO)注释[22],通过WEGO软件对所得到的unigenes做GO功能分类统计,从宏观上认识该物种的基因功能分布特征[23]

1.5 蛋白质编码框预测

将unigenes按照四大蛋白质数据库(NR、Swiss-Prot、COG和KEGG)的优先级顺序做Blast比对分析(E值<1×10-5),确定其unigenes编码区的核酸序列(序列方向5′→3′)和氨基酸序列,与四大蛋白质数据库未能比对上的unigenes,用ESTScan预测其蛋白质编码区(CDS)及氨基酸序列。

1.6 SSR与SNP分析

以组装出来的unigenes作为参考序列,使用MicroSAtellite(MISA)软件找出所有简单重复序列(Simple sequence repeat,SSR)并进行分析。将测序获得的转录组reads比对至参考序列,随后基于比对结果利用SOAPsnp检测样品中的单核苷酸多态性(Single nucleotide polymorphism,SNP)。

1.7 不同组织器官差异表达基因分析

为识别不同组织样本间表达水平存在显著差异的基因,使用DESeq对5种组织的unigenes进行两两比较[24],筛选标准为False discovery rate (FDR)≤0.001且|log2FoldChange|≥1(即表达量变化倍数≥2倍,上调或下调)。在获得差异表达基因后,进一步选取叶柄与根状茎两种组织的比较结果进行差异表达量、GO功能注释以及KEGG代谢通路富集分析。

2 结果与分析 2.1 转录组测序与组装

粗茎鳞毛蕨5种组织的转录组测序共获得242 779 056条raw reads,滤过后获得clean reads 233 883 608条,包含总碱基数35 082 541 200个。5种组织样本测序质控良好,clean reads质量合格,Q20(质量值≥20的碱基占比)均>96%,Q30(质量值≥30的碱基占比)均>92%,GC量均 < 50%(表 1)。

表 1 测序产量统计 Table 1 Statistics of sequencing yields
组织样品
Tissue sample
Raw reads条数
Number of raw reads
Clean reads条数
Number of clean reads
Clean reads碱基数/nt
Number of clean nucleotides/nt
质量值≥20碱基占比/%
Q20/%
质量值≥30碱基占比/%
Q30/%
GC含量/%
GC content/%
Leaf 47 701 822 45 537 714 6 830 657 100 97.50 92.20 49.78
Petiole1 51 831 394 49 554 358 7 433 153 700 97.14 93.85 49.77
Petiole2 45 878 090 44 765 978 6 714 896 700 97.01 93.04 47.94
Root 48 277 070 46 791 282 7 018 692 300 96.58 92.79 48.30
Stem 49 090 680 47 234 276 7 085 141 400 97.09 93.22 49.92

Trinity拼装得到的转录本去除冗余后,共获得147 690个unigenes,平均长度859 nt,N50为1 554 nt。长度为300 nt的unigenes数量最多,54 370个,显著高于其他长度区间的unigenes数量。随着序列长度的增加,unigenes的数量逐渐下降。在长序列区域,3 000 nt以上的unigenes数量明显减少,仅为5 199个。整体来看,短序列在unigenes数据集中占据主导地位,而较长的序列数量相对较少,尤其是大于2 000 nt的序列显著减少(图 1)。组装质量和各组织contig和unigenes长度分布统计详见补充材料。

图 1 粗茎鳞毛蕨unigenes长度分布 Fig. 1 Length distribution of D.crassirhizoma unigenes

2.2 Unigenes功能注释

在NR、NT、Swiss-Prot、KEGG、COG和GO数据库中分别注释到unigenes 91 758个(62.13%)、57 578个(38.99%)、60 535个(40.99%)、56 670个(38.37%)、41 446个(28.06%)、63 027个(42.68%)。

2.2.1 NR相似性分类

在NR数据库相似序列匹配度较高的物种中,葡萄(Vitis vinifera)所占比例最高,12 788条(13.94%);其次为桃(Prunus persica)11 614条(12.66%)、小立碗藓(Physcomitrella patens)6 937条(7.56%)、江南卷柏(Selaginella moellendorffii)6 741条(7.35%)、毛果杨(Populus tristis)6 553条(7.14%)、蓖麻(Ricinus communis)6 019条(6.56%)、巨云杉(Picea sitchensis)5 633条(6.14%)。其他物种比例均在6.14%以下,总占比为38.65%(图 2)。

图 2 粗茎鳞毛蕨转录组unigenes在NR数据库中的物种匹配分布 Fig. 2 Species distribution of transcriptome unigenes of D.crassirhizoma in NR database

2.2.2 COG功能分类

COG是对基因产物进行直系同源分类的数据库[25],将unigenes和COG数据库进行比对,可以宏观地认识粗茎鳞毛蕨整体的基因功能分布特征,预测unigenes可能的功能并对其做功能分类统计。结果显示共有41 446个unigenes被注释并归类到25个不同的功能类群中。其中,一般功能预测(R: general function prediction only)的基因数量最多,为12 011条;其次是转录(K: transcription)6 154条,复制、重组和修复(L: replication, recombination and repair)5 860条,其他种类基因丰度不尽相同,详见图 3

A: RNA processing and modification; B: chromatin structure and dynamics; C: energy production and conversion; D: cell cycle control, cell division, chromosome partitioning; E: amino acid transport and metabolism; F: nucleotide transport and metabolism; G: carbohydrate transport and metabolism; H: coenzyme transport and metabolism; I: lipid transport and metabolism; J: translation, ribosomal structure and biogenesis; K: transcription; L: replication, recombination and repair; M: cell wall/membrane/envelope biogenesis; N: cell motility; O: posttranslational modification, protein turnover, chaperones; P: inorganic ion transport and metabolism; Q: secondary metabolites biosynthesis, transport and catabolism; R: general function prediction only; S: function unknown; T: signal transduction mechanisms; U: intracellular trafficking, secretion, and vesicular transport; V: defense mechanisms; W: extracellular structures; X: nuclear structure; Y: cytoskeleton. 图 3 粗茎鳞毛蕨转录组的COG分布 Fig. 3 COG distribution of the D.crassirhizoma transcriptome

2.2.3 GO功能分类

GO是一种标准化的基因功能分类系统[26]。在GO功能数据库中比对粗茎鳞毛蕨unigenes,注释到的63 027条unigenes其功能被归为3大类共56个分支。其中,生物过程大类中的细胞过程涉及基因最多,共计41 350条,其次是代谢过程40 886条;细胞组成大类中的细胞、细胞部分和细胞器涉及基因分别有45 835、45 756和36 264条;分子功能大类中涉及的基因主要集中在连接体(30 130条)和催化过程(31 434条)(图 4)。

图 4 粗茎鳞毛蕨转录组的GO分类 Fig. 4 GO classification of the D.crassirhizoma transcriptome

2.2.4 KEGG代谢途径分类

KEGG是系统分析基因产物在细胞中的代谢途径以及这些基因产物功能的数据库[27]。以KEGG数据库作为参考,将转录组中的数据分成细胞过程、环境信息处理、遗传信息处理、代谢有机系统,共5大分支129条代谢通路(包括生化代谢通路、植物病原体互作、植物激素信号转导、苯丙烷类生物合成、RNA降解等)。根据基因注释度的大小进行排序,代谢途径Metabolic pathways: ko01100位列第一,共注释到unigenes 12 708个,占比22.42%。其次是次生代谢途径Biosynthesis of secondary metabolites: ko01110,共注释到unigenes 6 415个(11.32%)。排前10的代谢通路信息详见图 5

图 5 粗茎鳞毛蕨转录组的KEGG代谢途径分类 Fig. 5 KEGG metabolic pathway classification of the D.crassirhizoma transcriptome

对KEGG代谢途径进一步分析发现,共有3 985条unigenes参与包括苯丙烷类、黄酮、类黄酮等生物合成相关在内的23条次生代谢途径,详见图 6。其中,注释度最大的为苯丙烷类物质生物合成Phenylpropanoid biosynthesis:ko00940,共注释876条unigenes,占比1.55%。黄酮类相关通路注释度表现较大,如黄酮类生物合成途径Flavonoid biosynthesis: ko00941,共注释524条unigenes,占比0.92%;黄酮醇与黄酮生物合成途径Flavone and flavonol biosynthesis: ko00944,共注释238条unigenes,占比0.42%。注释度最小的为甜菜碱生物合成途径Betalain biosynthesis:ko00965,仅注释4条unigenes,占比0.01%。

图 6 粗茎鳞毛蕨转录组次生代谢相关unigenes的KEGG通路注释统计 Fig. 6 KEGG secondary metabolism pathway annotation statistics for D.crassirhizoma unigenes

2.2.5 预测编码蛋白框

对所有unigenes的CDS进行分析,通过Blast比对共获得CDS序列94 731个。从整体来看,大多数CDS序列的长度较短,尤其是长度在200-600 nt区间的序列数量显著较多,其中200 nt的序列数量最多。随着序列长度增加,CDS序列数量逐渐减少,长于2 000 nt的CDS序列数量明显下降(图 7)。

图 7 比对上的CDS序列长度分布 Fig. 7 Length distribution of aligned CDSs

预测的CDS序列长度以较短的序列为主,长序列的数量相对较少。其中200 nt长度的CDS序列数量最多,明显高于其他长度区间。总体上随着序列长度增加,CDS序列数量逐渐减少。长于1 300 nt的CDS序列区间,仅有极少数存在(图 8)。

图 8 ESTScan预测的CDS序列长度分布 Fig. 8 Length distribution of CDSs predicted by ESTScan

2.3 SSR特征分析

对粗茎鳞毛蕨转录组unigenes进行SSR分析,共得到28 079个SSRs。其中,二核苷酸SSRs序列数量最多,有19 690个,占70.12%,重复类型以AG/CT的数量最多。其次是三核苷酸SSR序列,数量为4 355个,占15.10%,重复类型以AAG/CTT数量最多(图 9)。

图 9 SSR长度分布统计 Fig. 9 Length distribution of SSRs

2.4 SNP情况分析

在粗茎鳞毛蕨根、叶、叶柄、根状茎、叶轴检测到的SNP数量分别为55 409、23 139、21 246、17 084、28 280个。根的SNP数量最多,Transition类型的SNP数量在根部最多,达到34 803个,其中A-G和C-T类型的SNP数量分别为17 397个和17 406个。根状茎的SNP数量最少,各类型的SNP数量在这5种组织中均最少(表 2)。

表 2 单核苷酸多态性数量分类统计表 Table 2 Statistics of different types of SNPs
单核苷酸多
态性类型SNP type

Root

Leaf
叶柄
Petiole
根状茎
Stem
叶轴
Leaf axes
Transition 34 803 14 535 12 945 10 766 17 094
A-G 17 397 7 310 6 462 5 364 8 635
C-T 17 406 7 225 6 483 5 402 8 459
Transversion 20 606 8 604 8 301 6 318 11 186
A-C 5 112 2 097 2 050 1 488 2 939
A-T 5 844 2 091 1 984 1 562 2 389
C-G 4 518 2 236 2 131 1 636 2 809
G-T 5 132 2 180 2 136 1 632 3 049
Total 55 409 23 139 21 246 17 084 28 280

2.5 差异表达基因分析

对粗茎鳞毛蕨5种组织进行两两比较,共发现10组存在显著差异表达的基因,其中,以根和叶柄显著差异基因数量最多(81 036个:上调32 486个、下调48 550个),根和叶轴次之(78 989个:上调47 496个、下调31 493个),叶和叶柄最少(15 671个:上调9 922个、下调5 749个)(图 10)。

图 10 差异表达基因统计 Fig. 10 Statistics of differentially expressed genes

进一步对根状茎与叶柄之间基因表达量差异进行分析,结果显示表达上调的基因有13 829个,表达下调的基因有24 250个(图 11)。

图 11 叶柄与根状茎差异表达火山图 Fig. 11 Volcano plot of differentially expressed genes between petioles and rhizomes

GO富集分析表明叶柄和根状茎之间差异基因共有18 967个,多数表现在细胞组成和分子功能上,其中以catalytic activity(催化活性:5 394个,28.44%)最多,其次为cytoplasm(细胞质:5 247个,27.66%)、cytoplasmic part(胞质部分:4 642个,24.47%)(图 12)。叶柄与根状茎富集到苯丙烷类生物合成途径(Phenylpropanoid biosynthesis)上的差异基因共有226个(2.13%),其中黄酮类合成基因占比最高(58个,25.66%),之后依次是木质素合成基因(52个,23.01%)、核心催化酶基因(41个,18.14%)、转录调控因子(27个,11.95%)、修饰酶基因(23个,10.18%)、转运蛋白基因(17个,7.52%),香豆素/木脂素合成基因数量最少(8个,3.54%)。其余代谢通路富集基因数量以代谢途径Metabolic pathways(2 623,24.68%)最多,详见图 13

图 12 叶柄和根状茎差异基因的GO功能分类 Fig. 12 GO functional classification of differentially expressed genes between petioles and rhizomes

图 13 叶柄和根状茎的代谢通路富集 Fig. 13 Metabolic pathway enrichment in petioles and rhizomes

3 讨论

高通量测序技术近些年发展迅速,已取得了重大进展,并被广泛应用于植物研究领域[28]。从测序质量结果来看,粗茎鳞毛蕨的转录组数据Q20均>96%,Q30>92%,质量整体表现良好,该结果符合已有研究的观点,即高质量测序数据通常要求Q20和Q30分别≥95%、≥90%[29]。良好的质量控制可以显著提升unigenes的覆盖度和准确性,粗茎鳞毛蕨转录组unigenes数量庞大,能满足转录组分析要求,且短序列在转录组中占据主导地位,长序列较为稀少。unigenes在各大数据库的注释结果显示,各数据的注释比例和已报道的植物[如黄三七[16]和独脚金(Striga asiatica)等[30]]类似。此外,粗茎鳞毛蕨转录组中存在大量序列特征及功能尚未知的unigenes,有可能是现有的数据库覆盖面不全,缺少相应的已知基因注释。

将unigenes序列与COG数据库进行比对,共获得25个功能类群,涵盖了多种与基本生命活动相关的生物学过程,从宏观层面揭示了粗茎鳞毛蕨基因的功能分布特征。GO分类揭示了粗茎鳞毛蕨根、根状茎在内的5种组织的转录组特性与生物过程、分子功能、细胞组成相关。通过对比KEGG数据库进一步研究粗茎鳞毛蕨基因在生物学上的功能,共得到5大分支129个KEGG标准代谢通路,包括生化代谢、植物激素信号转导等,这些可能与粗茎鳞毛蕨的呼吸作用、光合作用、水分和矿物质吸收等生理活动相关。此外,发现大量unigenes参与了苯丙烷类、黄酮和类黄酮等23条次生代谢相关的标准合成通路。酚类物质的合成在植物中主要与苯丙烷类代谢途径和莽草酸途径相关[31]。大量的unigenes参与了苯丙烷类和黄酮类代谢,这一结果不仅与粗茎鳞毛蕨主要成分为间苯三酚类的报道一致[8],也进一步验证了前人关于蕨类植物次生代谢产物主要由多酚类和黄酮类物质组成的研究结论[32-33]。这些发现为揭示粗茎鳞毛蕨间苯三酚类化合物的合成机制提供了有力的分子依据。通过粗茎鳞毛蕨的编码蛋白框分析得到的CDS序列,能够为后续新基因和蛋白质编码区的识别提供重要参考。

本研究共鉴定出28 079个SSRs位点,其中以二碱基重复类型最为丰富(19 690个,占比70.12%), 其次为三碱基重复类型(4 355个,占比15.51%)。在所有重复类型中,AG/CT和AAG/CCT的重复单元占比最高。该结果与越南安息香(Styrax tonkinensis)[34]、白芷(Angelica dahurica)[35]等植物的转录组SSR分析结果一致,进一步印证了前人关于大多数植物SSR重复类型以二碱基和三碱基重复为主的普遍观点,同时也反映不同植物中SSR类型分布仍存在一定差异性[36-38]。同一植物不同组织中的SNP分布可能反映基因组结构的变异,例如拷贝数变异和染色体结构变异。使用SOAPsnp检测粗茎鳞毛蕨5种组织的SNP,并对检测到的变异进行系统分类与统计分析,有助于深入理解基因组结构如何影响基因功能。

植物在生长过程中的生理、病理以及特定环境刺激的反应表现在功能基因的差异表达上[39],这些差异基因往往与该组织的发育、成熟和维持其功能密切相关。通过对粗茎鳞毛蕨5种组织差异表达基因两两对比分析,发现不同组织之间的差异基因数量相差较大。目前《中华人民共和国药典(一部)》(2020年版)[4]仅记载粗茎鳞毛蕨以干燥叶柄残基和根状茎入药,其叶柄和根状茎之间药效成分是否存在差异并没有明确指出。因此,本研究重点对叶柄和根状茎的差异基因进行分析,叶柄和根状茎之间存在较大差异,尤其是在某些基因的表达上,根状茎有更多的基因被显著下调,这可能与根状茎的功能、代谢需求及其在植物中的角色变化有关。叶柄和根状茎差异基因的GO富集结果表明,多数差异基因主要表现在细胞组成和分子功能上,类似情况也出现在其他植物中[40-41],这意味着这些基因在调节植物细胞的结构、代谢功能和应激反应过程中发挥着重要作用。粗茎鳞毛蕨叶柄和根状茎的代谢通路富集结果为后续进一步研究叶柄和根状茎之间药效成分差异提供了基础资料,并有助于挖掘药效成分合成的关键酶基因。

4 结论

本研究对粗茎鳞毛蕨5种组织进行转录组测序分析,共获得147 690个unigenes,其中在NR、NT、Swiss-Prot、KEGG、COG和GO数据库中分别有91 758、57 578、60 535、56 670、41 446和63 027个unigenes获得功能注释,系统揭示了粗茎鳞毛蕨的整体转录组组成与表达特征。对获得的unigenes进行GO分类统计,发现其广泛覆盖了分子功能、细胞组分和生物过程三大类的56个分支,同时在KEGG数据库中注释到129条标准通路,涵盖了23条次生代谢通路。共鉴定预测到94 731个蛋白编码框序列、28 079个SSRs位点以及不同组织间数量丰富的SNP变异。差异表达分析显示,叶柄与根状茎间存在大量表达变化(上调基因13 829个、下调基因24 250个),其中富集到226个与苯丙烷类代谢途径相关的差异基因。上述研究结果为后续进一步研究粗茎鳞毛蕨次生代谢通路、间苯三酚类、黄酮类等物质的生物合成、基因功能鉴定等,提供了转录组相关的数据基础,同时促进了粗茎鳞毛蕨资源的合理应用。

参考文献
[1]
中国科学院中国植物志编辑委员会. 中国植物志: 第五卷第一分册[M]. 北京: 科学出版社, 2000: 149.
[2]
(清)孙星衍辑注. 徐斌校注. 神农本草经精注易读本[M]. 北京: 中国中医药出版社, 2019: 198.
[3]
楼之岑, 秦波. 常用中药材品种整理和质量研究北方编第1册[M]. 北京: 北京医科大学、中国协和医科大学联合出版社, 1995: 57.
[4]
国家药典委员会. 中华人民共和国药典: 一部[S]. 2020版. 北京: 中国医药科技出版社, 2025: 353-354.
[5]
王研加. 粗茎鳞毛蕨化学成分鉴定及其生物活性研究[D]. 哈尔滨: 哈尔滨师范大学, 2023.
[6]
魏凯欣, 宋雄辉, 刘向前, 等. 绵马贯众化学成分及其抗炎活性研究[J]. 天然产物研究与开发, 2025, 37(1): 65-73.
[7]
WANG J, YAN Y T, FU S Z, et al. Anti-influenza virus (H5N1) activity screening on the phloroglucinols from rhizomes of Dryopteris crassirhizoma[J]. Molecules, 2017, 22(3): 431. DOI:10.3390/molecules22030431
[8]
YUK H J, KIM J Y, SUNG Y Y, et al. Phloroglucinol derivatives from Dryopteris crassirhizoma as potent xanthine oxidase inhibitors[J]. Molecules, 2021, 26(1): 122.
[9]
吴以岭. 解读连花清瘟胶囊[J]. 中国医药指南, 2005, 3(11): 120-121.
[10]
王光月, 姜璐璐, 霍若琳, 等. 连花清瘟胶囊定性鉴别优化研究[J]. 吉林中医药, 2021, 41(12): 1669-1672.
[11]
PIÉTU G, MARIAGE-SAMSON R, FAYEIN N A, et al. The Genexpress IMAGE knowledge base of the human brain transcriptome: a prototype integrated resource for functional and computational genomics[J]. Genome Research, 1999, 9(2): 195-209. DOI:10.1101/gr.9.2.195
[12]
陈士林, 朱孝轩, 李春芳, 等. 中药基因组学与合成生物学[J]. 药学学报, 2012, 47(8): 1070-1078.
[13]
张晓萌, 李健春, 王琼, 等. 转录组测序技术在中医药领域的应用[J]. 中国现代中药, 2016, 18(8): 1084-1087.
[14]
王思冉, 李慧, 阿日查, 等. 杜仲叶片不同发育时期转录组及绿原酸相关基因分析[J]. 中南林业科技大学学报, 2025, 45(7): 174-187.
[15]
胡菊, 杨秀玲, 张恒, 等. 朱砂根及其变种红凉伞叶片转录组测序特征分析[J]. 西南林业大学学报(自然科学), 2025, 45(6): 108-117.
[16]
李依民, 彭亮, 杨冰月, 等. 基于高通量测序技术的黄三七根茎转录组数据分析[J]. 中草药, 2018, 49(21): 4983-4990.
[17]
陈华园, 庞伟灿, 李翠, 等. 岩黄连全长转录组测序及生物信息学分析[J]. 种子, 2024, 43(12): 27-33, 39.
[18]
GRABHERR M G, HAAS B J, YASSOUR M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome[J]. Nature Biotechnology, 2011, 29(7): 644-652. DOI:10.1038/nbt.1883
[19]
BOECKMANN B, BAIROCH A, APWEILER R, et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003[J]. Nucleic Acids Research, 2003, 31(1): 365-370. DOI:10.1093/nar/gkg095
[20]
GÖTZ S, GARCÍA-GÓMEZ J M, TEROL J, et al. High-throughput functional annotation and data mining with the Blast2GO suite[J]. Nucleic Acids Research, 2008, 36(10): 3420-3435. DOI:10.1093/nar/gkn176
[21]
PRUITT K D, TATUSOVA T, MAGLOTT D R. NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins[J]. Nucleic Acids Research, 2005, 33: D501-D504. DOI:10.1093/nar/gki476
[22]
CONESA A, GÖTZS S, GARCÍA-GÓMEZ J M, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research[J]. Bioinformatics, 2005, 21(18): 3674-3676. DOI:10.1093/bioinformatics/bti610
[23]
ASHBURNER M, BALL C A, BLAKE J A, et al. Gene ontology: tool for the unification of biology[J]. Nature Genetics, 2000, 25(1): 25-29. DOI:10.1038/75556
[24]
LOVE M I, HUBER W, ANDERS S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J]. Genome Biology, 2014, 15(12): 550. DOI:10.1186/s13059-014-0550-8
[25]
TATUSOV R L, FEDOROVA N D, JACKSON J D, et al. The COG database: an updated version includes eukaryotes[J]. BMC Bioinformatics, 2003, 4: 41. DOI:10.1186/1471-2105-4-41
[26]
LIU X B, MA L, ZHANG A H, et al. High-throughput analysis and characterization of Astragalus membranaceus transcriptome using 454 GS FLX[J]. PLoS One, 2014, 9(5): e95831. DOI:10.1371/journal.pone.0095831
[27]
KANEHISA M, GOTO S, KAWASHIMA S, et al. The KEGG resource for deciphering the genome[J]. Nucleic Acids Research, 2004, 32(suppl_1): D277-D280.
[28]
HAILE Z M, NAGPALA-DE GUZMAN E G, MORETTO M, et al. Transcriptome profiles of strawberry (Fragaria vesca) fruit interacting with Botrytis cinerea at different ripening stages[J]. Frontiers in Plant Science, 2019, 10: 1131. DOI:10.3389/fpls.2019.01131
[29]
WATSON M. Quality assessment and control of high-throughput sequencing data[J]. Frontiers in Genetics, 2014, 5: 235.
[30]
张龙, 蔡时可, 孙金锴, 等. 独脚金转录组特征分析[J]. 广州中医药大学学报, 2023, 40(8): 2042-2048.
[31]
VOGT T. Phenylpropanoid biosynthesis[J]. Molecular Plant, 2010, 3(1): 2-20. DOI:10.1093/mp/ssp106
[32]
BANDYOPADHYAY A, DEY A. Medicinal pterido-phytes: ethnopharmacological, phytochemical, and clinical attributes[J]. Beni-Suef University Journal of Basic and Applied Sciences, 2022, 11(1): 113. DOI:10.1186/s43088-022-00283-3
[33]
IBADULLAYEVA S J, MOVSUMOVA N V, SHIRA-LIYEVA G S, et al. Pteridophytes: ethnobotanical use and active chemical composition[J]. Indian Journal of Traditional Knowledge, 2022, 21(2): 353-359.
[34]
王小敏, 许婳婳, 詹若挺, 等. 基于转录组测序的越南安息香根、茎和叶基因表达分析[J]. 中草药, 2021, 52(8): 2392-2399.
[35]
CHEN C, CHEN Y J, HUANG W J, et al. Mining of simple sequence repeats (SSRs) loci and development of novel transferability-across EST-SSR markers from de novo transcriptome assembly of Angelica dahurica[J]. PLoS One, 2019, 14(8): e0221040. DOI:10.1371/journal.pone.0221040
[36]
AN Y L, XIA X B, ZHENG H Y, et al. Multi-genome comprehensive identification of SSR/SV and development of molecular markers database to serve Sorghum bicolor (L.) breeding[J]. BMC Genomic Data, 2023, 24(1): 62. DOI:10.1186/s12863-023-01165-y
[37]
TAO A F, LI Y Q, CHEN J H, et al. Development of Roselle (Hibiscus sabdariffa L.) transcriptome-based simple sequence repeat markers and their application in Roselle[J]. Plants, 2024, 13(24): 3517. DOI:10.3390/plants13243517
[38]
PARK H, HEO T H, CHO J, et al. Evaluation and characteristic analysis of SSRs from the transcriptomic sequences of Perilla crop (Perilla frutescens L.)[J]. Gene, 2025, 933: 148938. DOI:10.1016/j.gene.2024.148938
[39]
FANG Y, HAN Y P, FANG Y J, et al. Epigenetic regulation modulates seasonal temperature-dependent growth of soybean in southern China[J]. Plant Biotechnology Journal, 2025, 23(10): 4580-4601. DOI:10.1111/pbi.70243
[40]
HUANG L B, ZHANG L Y, ZHANG P, et al. Comparative transcriptomes and WGCNA reveal hub genes for spike germination in different quinoa lines[J]. BMC Genomics, 2024, 25(1): 1231. DOI:10.1186/s12864-024-11151-y
[41]
BAI X D, ZHENG Y, CAO L, et al. Transcriptome and metabolome conjoint analysis revealed that PaGLK affects photosynthesis and composition of root exudates in poplar[J]. Plant Molecular Biology Reporter, 2025, 43(3): 1369-1379. DOI:10.1007/s11105-025-01550-0