幼虫、青年和成年日本医蛭不同生长阶段差异表达基因比较分析

2022-02-03闫秀营任可马晓聪单会荃郑雅琳刘庆友周维官郑景辉

南方农业学报 2022年10期

闫秀营，任可，马晓聪，单会荃，郑雅琳，刘庆友*，周维官，郑景辉*

（1广西大学动物科学技术学院/亚热带农业生物资源保护与利用国家重点实验室，广西南宁 530004；2广西中医药大学附属瑞康医院，广西南宁 530011；3泰国天然水蛭素股份有限公司，泰国曼谷 10700）

0 引言

【研究意义】日本医蛭作为传统中药已有几千历史，是《中国药典》收录的唯一一个吸食血液的水蛭（王懿等，2021），含有最强凝血酶天然抑制剂——水蛭素（Koeppen et al.，2014），临床上可用于治疗多种疾病，包括高血压、高脂血、脑梗死及恶性肿瘤等（黄慧等，2021；李晓文等，2021；王懿等，2021）。除此之外，水蛭的免疫系统含有多种具有杀菌作用的肽，可用于研制新一代的抗菌药物（余米等，2021b）。水蛭的身体和中枢神经系统均能再生，且其神经细胞结构和作用与人类的神经相似，因此常被用作神经生物学领域研究的模式生物（Meiser et al.，2019；Cohen et al.，2020；Heath-Heckman et al.，2021）。近年来，随着日本医蛭生存环境的破坏和捕猎数量的增加，野生水蛭数量逐渐减少。日本医蛭个体小、繁殖慢，养殖经济效益远不如其他水蛭（杨谋等，2018；余米等，2019，2021a），其数量已无法满足市场需求。因此，开展日本医蛭转录组分析，了解其遗传信息及生长发育的调控机制，对日本医蛭基因资源的开发利用及药用品质提升具有重要意义。【前人研究进展】日本医蛭从幼苗发育生长到成虫体型大约需要80 d（张海生，2020），但成熟周期较长、繁殖慢，3年生日本医蛭才开始产卵。与宽体金线蛭和菲牛蛭相比，日本医蛭的平均产卵茧数最少（马燕，2016）。美国食品和药物管理局（Food and Drug Administration，FDA）于2004年6月批准医用水蛭可作为医疗设备予以应用（Singh and Rajoria，2020），具体是指水蛭疗法，即将水蛭含有许多生物活性物质的唾液注入生物体内（Koeppen et al.，2020），用于皮瓣改善、伤口愈合、重建手术、急性静脉充血恢复，以及对脑缺血再灌注损伤的屏蔽作用（Mumcuoglu，2014；Dong et al.，2016；Mousavi et al.，2016；Sig et al.，2017）。水蛭唾液中最具代表性的水蛭素是一种酸性单链多肽，由64~66个氨基酸残基组成，含有多个二硫键，能与凝血酶结合，从而抑制其与纤维蛋白原的结合，发挥抗凝和溶栓作用（王斌等，2017；杨谋等，2018；侯宁，2021）。随着高通量测序技术的发展，对水蛭分子生物学方面的研究也越来越深入全面。Liu等（2018）对3种不同种水蛭（棒纹牛蛭、宽体金线蛭和洞穴山蛭）进行比较转录组学研究，结果发现3种水蛭的嗅觉转导通路基因存在明显差异；邢月婷等（2020）对饥饿状态和取食的日本医蛭进行转录组分析，结果发现许多抗凝和抗炎症物质；余米等（2021a）基于高通量测序的日本医蛭转录组分析结果表明，不同温度处理下的日本医蛭转录组和简单重复序列标记变化较快，且水蛭素、FBXO家族和Dmrt等抗凝血及繁殖相关基因变化较明显；Tong等（2022）通过比较宽体金线蛭、日本医蛭和淡水水蛭的基因组，结果发现3种水蛭在同源基因、基因共线性及基因家族方面存在明显差异。【本研究切入点】至今，有关日本医蛭不同生长阶段的组学研究鲜见报道，进而制约其基因资源的高效开发利用。【拟解决的关键问题】基于RNA-Seq技术对不同生长发育阶段的日本医蛭进行转录组测序分析，并对所得基因进行功能注释和分类，旨在筛选出大量差异表达基因（Differentially expressed genes，DEGs），为后期开展日本医蛭养殖、生长发育及相关功能基因研究提供科学依据。

1 材料与方法

1.1 样品采集

供试用日本医蛭均采集自钦州市展宏中药材公司，样品分为幼虫期（孵化10 d以内，HNL）、青年虫期（孵化后2个月左右，HNY）和成年虫期（孵化后2年左右，HNA）。根据日本医蛭生理学解剖结构将其分为幼虫前后段、青年虫前后段、成虫口吸盘、性腺、尾吸盘和剩余虫体等24个样本（表1），经液氮迅速冷冻后置于-80℃冰箱保存备用。同时，另取9个样品进行组织学观察，固定于4%多聚甲醛中，4℃保存，每组取3个个体作为生物学重复。

表1 样品编号Table 1 Sample number

1.2 组织形态学观察

经常规石蜡包埋、切片及苏木精—伊红染色后，通过显微镜（EVOS FL Auto）观察日本医蛭的组织学特征，并以显微摄影系统（EVOS Auto）采集图像。

1.3 cDNA文库构建及序列数据分析

委托安诺优达基因科技（北京）有限公司进行样品RNA提取及完成cDNA文库构建，通过Illumina HiSeq 4000平台进行测序，获取长度为150 bp的双末端序列。使用Fastp对每个样品进行质量控制（Chen et al.，2018），通过HISAT2将各样本的Reads映射到参考基因组，以SAMtools（v1.11）对SAM文件进行排序并转换为BAM文件（Li et al.，2009）。然后应用StringTie进行组装和去冗余（Pertea et al.，2015），获得最后的转录本；运用GffCompare对转录本与GTF中的基因注释信息进行比对（Pertea and Pertea，2020），以StringTie的选项“-e-B-P20”估计转录本丰度，丰度结果以“.balltown”结尾的文件呈现，并使用prepDE.py比较这些文件。

1.4 差异表达基因筛选

采用FPKM法计算各样本中所有基因的表达量，通过DESeq2筛选不同组间的差异表达基因。筛选原则：错误发现率（FDR）＜0.05和|log2FC|＞1的基因即被认为是不同生长阶段间的差异表达基因。

1.5 样本时间序列分析

采用TCseq分析3个生长阶段的差异表达基因并进行聚类分析，直观显示不同生长阶段的不同基因表达模式（Ma et al.，2021）。

1.6 差异表达基因功能富集分析

通过clusterProfiler分别进行差异表达基因的GO功能注释分析和KEGG信号通路富集分析（Yu et al.，2012），FDR＜0.05即为显著富集。

2.1 不同生长阶段日本医蛭的发育性状及组织学观察结果

经4%多聚甲醛固定的日本医蛭幼虫（L）、青年虫（Y）和成虫（A），如图1-A所示。计算不同生长阶段日本医蛭的体质量，其中，幼虫体质量为0.136 g、青年虫体质量为1.184 g、成年虫体质量为1.524 g，不同生长阶段日本医蛭的体质量存在极显著差异（P＜0.01）（图1-B）。此外，不同生长阶段日本医蛭的石蜡切片观察结果（图1-C）显示，日本医蛭成年虫期切片中的肌纤维明显多于幼虫期和青年虫期。

图1 日本医蛭的表型数据及石蜡切片观察结果Fig.1 Phenotypic data and paraffin section of H.nipponia

2.2 RNA测序分析结果

构建24个日本医蛭cDNA文库，将不同生长阶段不同部位的数据进行组合，从日本医蛭幼虫期、青年虫期和成年虫期各获得3组样品（HNL1~HNL3、HNY1~HNY3和HNA1~HNA3），共产生541124257条原始序列（Raw reads）；经Fastp质量控制后，获得517805628条有效序列（Clean reads）。所有样本的GC含量十分接近，位于41%~45%之间（表2）。以上结果表明测序数据质量较好，可满足后续分析的要求。

表2 日本医蛭转录组测序数据统计结果Table 2 Summary statistics of H.nipponia transcriptome sequencing

RNA-Seq共组装出20430个mRNAs转录本，其中88个为新发现的mRNAs转录本。不同样本的基因表达情况如图2-A所示。主成分分析（PCA）结果表明，9个样品的基因表达量数据分布基本一致（图2-B），说明不同分组样本间可较清晰地划分，组内相关性也较高，可为后续分析提供可靠数据。

图2 mRNAs表达测序分析结果Fig.2 mRNAs expression sequencing analysis

2.3 差异表达基因及样本时间序列分析结果

根据FDR＜0.05和|log2FC|＞1的差异表达基因筛选标准，对日本医蛭不同生长阶段进行两两对比，结果共获得1348个差异表达基因。其中，在幼虫期—青年虫期（HNL vs HNY）中获得238个差异表达基因，在幼虫期—成虫期（HNL vs HNA）中获得976个差异表达基因，在青年虫期—成虫期（HNY vs HNA）中获得763个差异表达基因（图3）。可见，成虫期较其他2个生长期的基因表达差异明显，表现活跃；而幼虫和青年虫间的差异相对较小，表现较稳定，说明生长发育阶段对日本医蛭基因表达影响显著。进一步分析差异表达基因的表达趋势，结果（图4）发现基因表达趋势表现为：3个样本（Cluster 2、Cluster 4和Cluster 6）持续下降、2个样本（Cluster 3和Cluster 5）持续上升、3个样本（Cluster 1、Cluster 7和Cluster 8）青年期高表达和1个样本（Cluster 9）青年期低表达。

图3 日本医蛭差异表达基因统计结果Fig.3 Statistics of H.nipponia DEGs

图4 日本医蛭差异表达基因样本时间序列分析结果Fig.4 DEGs sample time series analysis of H.nipponia

2.4 差异表达基因样本时间序列GO功能注释和KEGG信号通路富集分析结果

为揭示调控日本医蛭发育的分子机制，对其上调趋势的Cluster 4（含329个差异表达基因）和下调趋势的Cluster 5（含176个差异表达基因）进行GO功能注释分析，结果（图5）表明，Cluster 4主要注释在肌球蛋白复合体、运动活动、肌动蛋白丝结合、蛋白运输及金属酶等功能条目上，Cluster 5主要注释于离子转运、信号传导及氧运输等功能条目上。KEGG信号通路富集分析结果（图6）显示，差异表达基因样本时间序列主要富集在氨基酸代谢、苯丙氨酸、酪氨酸和色氨酸生物合成及原核生物碳固定等通路上，但富集结果显著性不高（FDR＞0.05）。

图5 差异表达基因样本时间序列GO功能注释分析结果Fig.5 GO functional annotation analysis of sample time series of DEGs

图6 差异表达基因样本时间序列KEGG信号通路富集分析结果Fig.6 KEGG signal pathway enrichment analysis of sample time series of DEGs

2.5 不同生长阶段差异表达基因GO功能注释和KEGG信号通路富集分析结果

对不同生长阶段差异表达基因进行GO功能注释分析和KEGG信号通路富集分析，结果发现，幼虫期与成虫期相比，上调差异表达基因显著注释在肌球蛋白复合体、细胞骨架、金属肽酶活性和金属内肽酶活性等GO功能条目上（图7-A），下调差异表达基因主要注释在几丁质代谢过程、对伤害的反应、氧结合、血红素结合和离子运输等GO功能条目上（图7-B）；上调差异表达基因主要富集于氨基酸代谢、糖胺聚糖生物合成—硫酸软骨素/硫酸皮肤素等KEGG信号通路上（图7-C），但富集不显著（FDR＞0.05），下调差异表达基因则显著富集在丙烯酰胺转移酶、卟啉与叶绿素代谢、苯丙氨酸/酪氨酸/色氨酸生物合成等KEGG信号通路上（图7-D）。此外，其他不同生长阶段的差异表达基因的GO功能注释和KEGG信号通路富集分析结果显示，上调的代谢通路主要富集在能量代谢和肌肉发育方面，而下调的通路富集在代谢和信号传导方面（表3）。这可能是日本医蛭幼虫在生长过程中大量摄食并通过消化、代谢和合成物质为自身发育提供能量，而成年日本医蛭发育成熟后其脱皮现象已完成，与其生长发育相关的几丁质代谢过程、生物体必需氨基酸合成通路均下调，尤其是成年日本医蛭已适应环境，对于环境刺激的敏感度也有所降低。

表3 差异表达基因显著富集的KEGG信号通路Table 3 KEGG signaling pathway with significant enrichment of DEGs

图7 日本医蛭幼虫期—成年期差异表达基因的功能富集分析结果Fig.7 The DEGs enrichment analysis in HNL vs HNA periods of H.nipponia

2.6 日本医蛭发育性状相关的潜在基因

进一步对差异表达基因的富集途径分析发现，GO功能注释获得肌球蛋白复合体、细胞骨架等组分中的差异表达基因在日本医蛭的生长发育过程中（HNL→HNY→HNA）呈上调趋势（图8-A）；而表面蛋白（SPA17）、IQ motif containing G（IQCG）、ACTB、5-羟色胺（5-HT）和HOX基因家族等同源基因可能参与日本医蛭的性腺发育或产卵调控（图8-B）；参与神经发育相关的微管蛋白、细丝蛋白A（FLNA）、原钙黏蛋白（PCDHGA10）、溶质载体家族5成员8（SLC5A8）和NK型同源盒基因转录因子2.8（NKX2-8）等基因在各生长阶段均有表达，且以在幼虫期的表达量最高（图8-C）。本研究还发现，与抗凝相关的蛋白基因在青年虫中的表达量较高，包括水蛭素（Hirudin）、抗蛋白酶（ANTA）、金属蛋白酶家族（ADAMTS）、Therostasin和血小板糖蛋白抑制剂（DECO）等（图8-D）。

图8 日本医蛭生长发育相关差异表达基因的表达情况Fig.8 The expression of DEGs related to growth and development of H.nipponia

3 讨论

本研究对不同所长阶段的日本医蛭进行转录组学分析，以期为研究日本医蛭发育过程中的生理机制及提高其人工养殖效益提供新的理论依据。在3个生长阶段中共发现1348个差异表达基因，且差异表达基因数量随着日本医蛭的发育而逐渐减少，可能与3月龄青年虫已发育至成虫有关。差异表达基因的GO功能注释分析及KEGG信号通路富集分析结果表明，肌球蛋白复合体、细胞骨架和细胞生长发育相关信号通路呈上调趋势，是由于幼虫和青年虫时期日本医蛭生长速度快，需大量摄食以获得能量，与非生殖阶段（摄食阶段）成虫期相比观察到更积极摄食行为的结论（Xu et al.，2021）一致；下调差异表达基因则主要与信号转导、苯丙氨酸、酪氨酸和色氨酸生物合成等信号通路相关，究其原因是随着日本医蛭的发育，其对食物的需求度及代谢水平均逐渐下降。

已有研究表明，动物的比生长速率在不同生长阶段存在显著变化，说明潜在的肌肉生长调节机制可能不同（Zhao et al.，2021）。本研究发现，成年日本医蛭的肌纤维数量明显多于幼虫和青年虫，GO功能注释分析也显示HNL vs HNA和HNY vs HNA的上调差异表达基因主要富集在肌球蛋白复合体、细胞骨架等条目上。其中，副肌球蛋白是环节动物的肌肉组成蛋白，在肌肉紧张收缩机制中发挥重要作用（Sonobe et al.，2016；Li et al.，2021）。原肌球蛋白是肌肉细肌丝中与肌动蛋白结合的蛋白，除了完成基本的收缩功能外，还参与完成体内多种生命活动，如胞质分裂及信号转导等；通过与肌动蛋白的细肌丝螺旋形结合存在，而调节与肌球蛋白间的相互作用（Pavadai et al.，2020）。ACTB是肌动蛋白家族的重要成员之一，几乎参与所有形式的细胞和细胞器运动，在生长发育、生物调控及修复过程中发挥重要作用（Hammer et al.，2019）。肌球蛋白由MYL和MYH组成，其中，MYH家族是调节肌肉发育的关键成员之一，主要通过调节细胞增殖、分化来影响肌肉的生成（Liu et al.，2020；Hill et al.，2021），这些基因均在成年日本医蛭中高表达。

抗凝血能力是水蛭研究的热点之一。本研究发现，在日本医蛭体内除了水蛭素以外，抗凝血基因ANTA和DECO也有表达。水蛭素可抑制胰蛋白酶、糜蛋白酶、组织蛋白酶G和组织激肽释放酶，但不抑制凝血因子Xa活性（王斌等，2017；Lu et al.，2018）。相反，ANTA是凝血因子Xa的有效抑制剂。DECO能抑制纤维蛋白原表达，与糖蛋白IIb-IIIa复合物上的血小板受体相互作用，防止血液凝结（Kwak et al.，2019）。抗凝基因在日本医蛭青年期的表达量高于成年，该结论为青年日本医蛭的药用提供了理论依据。

水蛭是典型的雌雄同体，雄性腺优先于雌性腺发育成熟。与其他品种水蛭相比，日本医蛭的繁殖性能较低（马燕，2016）。本研究通过转录组测序分析不同生长发育阶段的差异表达基因，筛选出与性腺发育或卵发生相关的候选基因或蛋白。其中，SP17是一类与精子发生和形成、精卵识别及顶体反应等受精活动密切相关的蛋白（Minhas et al.，2020）。余米等（2021a）研究证实，在18℃下SP基因在日本水蛭体内高表达。IQCG也参与精子的发生过程，是纤毛/鞭毛运动的关键调节器，与精子运动有关（Harris et al.，2014；Pan et al.，2016）。关于卵巢发育的研究表明，5-HT在性腺成熟和连续产卵方面发挥多种生物学作用，是通过G蛋白偶联受体超家族5-HT受体发挥作用，调节无脊椎动物和脊椎动物的卵母细胞发育与成熟（Song et al.，2016）。HOX基因是一个调节分子家族且编码高度保守的转录因子，其表达受性类固醇的调控，在生殖道发育、子宫内膜周期生长和胚胎植入中发挥重要作用（Du and Taylor，2015；Kondo et al.，2017）。这些基因在日本医蛭成年期显著上调，可能与其繁殖性能相关。

日本医蛭是研究神经再生机制的理想模型系统，具有在整个生命周期内再生中枢神经系统的能力。此外，日本医蛭中枢神经元持续扩张，即其生长和增加突触耦合的机制可能永远不会被完全关闭（Meriaux et al.，2011；Cohen et al.，2020）。本研究发现，微管蛋白（刁磊等，2019）、FLNA（周翔宇，2021）及轴突生长诱向因子（Netrin）（孟凡伟等，2003）等与神经发育相关的基因在3个生长阶段均有表达，其中，FLNA、PCDHGA10（Schoch et al.，2017）和NKX2-8在幼虫期高表达，故可作为筛选日本医蛭的候选分子遗传标记。