閆秀營，任可，馬曉聰，單會荃，鄭雅琳，劉慶友*，周維官，鄭景輝*

（1廣西大學動物科學技術(shù)學院/亞熱帶農(nóng)業(yè)生物資源保護與利用國家重點實驗室，廣西南寧 530004；2廣西中醫(yī)藥大學附屬瑞康醫(yī)院，廣西南寧 530011；3泰國天然水蛭素股份有限公司，泰國曼谷 10700）

0 引言

【研究意義】日本醫(yī)蛭作為傳統(tǒng)中藥已有幾千歷史，是《中國藥典》收錄的唯一一個吸食血液的水蛭（王懿等，2021），含有最強凝血酶天然抑制劑——水蛭素（Koeppen et al.，2014），臨床上可用于治療多種疾病，包括高血壓、高脂血、腦梗死及惡性腫瘤等（黃慧等，2021；李曉文等，2021；王懿等，2021）。除此之外，水蛭的免疫系統(tǒng)含有多種具有殺菌作用的肽，可用于研制新一代的抗菌藥物（余米等，2021b）。水蛭的身體和中樞神經(jīng)系統(tǒng)均能再生，且其神經(jīng)細胞結(jié)構(gòu)和作用與人類的神經(jīng)相似，因此常被用作神經(jīng)生物學領(lǐng)域研究的模式生物（Meiser et al.，2019；Cohen et al.，2020；Heath-Heckman et al.，2021）。近年來，隨著日本醫(yī)蛭生存環(huán)境的破壞和捕獵數(shù)量的增加，野生水蛭數(shù)量逐漸減少。日本醫(yī)蛭個體小、繁殖慢，養(yǎng)殖經(jīng)濟效益遠不如其他水蛭（楊謀等，2018；余米等，2019，2021a），其數(shù)量已無法滿足市場需求。因此，開展日本醫(yī)蛭轉(zhuǎn)錄組分析，了解其遺傳信息及生長發(fā)育的調(diào)控機制，對日本醫(yī)蛭基因資源的開發(fā)利用及藥用品質(zhì)提升具有重要意義?！厩叭搜芯窟M展】日本醫(yī)蛭從幼苗發(fā)育生長到成蟲體型大約需要80 d（張海生，2020），但成熟周期較長、繁殖慢，3年生日本醫(yī)蛭才開始產(chǎn)卵。與寬體金線蛭和菲牛蛭相比，日本醫(yī)蛭的平均產(chǎn)卵繭數(shù)最少（馬燕，2016）。美國食品和藥物管理局（Food and Drug Administration，F(xiàn)DA）于2004年6月批準醫(yī)用水蛭可作為醫(yī)療設(shè)備予以應(yīng)用（Singh and Rajoria，2020），具體是指水蛭療法，即將水蛭含有許多生物活性物質(zhì)的唾液注入生物體內(nèi)（Koeppen et al.，2020），用于皮瓣改善、傷口愈合、重建手術(shù)、急性靜脈充血恢復(fù)，以及對腦缺血再灌注損傷的屏蔽作用（Mumcuoglu，2014；Dong et al.，2016；Mousavi et al.，2016；Sig et al.，2017）。水蛭唾液中最具代表性的水蛭素是一種酸性單鏈多肽，由64~66個氨基酸殘基組成，含有多個二硫鍵，能與凝血酶結(jié)合，從而抑制其與纖維蛋白原的結(jié)合，發(fā)揮抗凝和溶栓作用（王斌等，2017；楊謀等，2018；侯寧，2021）。隨著高通量測序技術(shù)的發(fā)展，對水蛭分子生物學方面的研究也越來越深入全面。Liu等（2018）對3種不同種水蛭（棒紋牛蛭、寬體金線蛭和洞穴山蛭）進行比較轉(zhuǎn)錄組學研究，結(jié)果發(fā)現(xiàn)3種水蛭的嗅覺轉(zhuǎn)導(dǎo)通路基因存在明顯差異；邢月婷等（2020）對饑餓狀態(tài)和取食的日本醫(yī)蛭進行轉(zhuǎn)錄組分析，結(jié)果發(fā)現(xiàn)許多抗凝和抗炎癥物質(zhì)；余米等（2021a）基于高通量測序的日本醫(yī)蛭轉(zhuǎn)錄組分析結(jié)果表明，不同溫度處理下的日本醫(yī)蛭轉(zhuǎn)錄組和簡單重復(fù)序列標記變化較快，且水蛭素、FBXO家族和Dmrt等抗凝血及繁殖相關(guān)基因變化較明顯；Tong等（2022）通過比較寬體金線蛭、日本醫(yī)蛭和淡水水蛭的基因組，結(jié)果發(fā)現(xiàn)3種水蛭在同源基因、基因共線性及基因家族方面存在明顯差異?！颈狙芯壳腥朦c】至今，有關(guān)日本醫(yī)蛭不同生長階段的組學研究鮮見報道，進而制約其基因資源的高效開發(fā)利用?！緮M解決的關(guān)鍵問題】基于RNA-Seq技術(shù)對不同生長發(fā)育階段的日本醫(yī)蛭進行轉(zhuǎn)錄組測序分析，并對所得基因進行功能注釋和分類，旨在篩選出大量差異表達基因（Differentially expressed genes，DEGs），為后期開展日本醫(yī)蛭養(yǎng)殖、生長發(fā)育及相關(guān)功能基因研究提供科學依據(jù)。

1 材料與方法

1.1 樣品采集

供試用日本醫(yī)蛭均采集自欽州市展宏中藥材公司，樣品分為幼蟲期（孵化10 d以內(nèi)，HNL）、青年蟲期（孵化后2個月左右，HNY）和成年蟲期（孵化后2年左右，HNA）。根據(jù)日本醫(yī)蛭生理學解剖結(jié)構(gòu)將其分為幼蟲前后段、青年蟲前后段、成蟲口吸盤、性腺、尾吸盤和剩余蟲體等24個樣本（表1），經(jīng)液氮迅速冷凍后置于-80℃冰箱保存?zhèn)溆谩Ｍ瑫r，另取9個樣品進行組織學觀察，固定于4%多聚甲醛中，4℃保存，每組取3個個體作為生物學重復(fù)。

表1 樣品編號Table 1 Sample number

1.2 組織形態(tài)學觀察

經(jīng)常規(guī)石蠟包埋、切片及蘇木精—伊紅染色后，通過顯微鏡（EVOS FL Auto）觀察日本醫(yī)蛭的組織學特征，并以顯微攝影系統(tǒng)（EVOS Auto）采集圖像。

1.3 cDNA文庫構(gòu)建及序列數(shù)據(jù)分析

委托安諾優(yōu)達基因科技（北京）有限公司進行樣品RNA提取及完成cDNA文庫構(gòu)建，通過Illumina HiSeq 4000平臺進行測序，獲取長度為150 bp的雙末端序列。使用Fastp對每個樣品進行質(zhì)量控制（Chen et al.，2018），通過HISAT2將各樣本的Reads映射到參考基因組，以SAMtools（v1.11）對SAM文件進行排序并轉(zhuǎn)換為BAM文件（Li et al.，2009）。然后應(yīng)用StringTie進行組裝和去冗余（Pertea et al.，2015），獲得最后的轉(zhuǎn)錄本；運用GffCompare對轉(zhuǎn)錄本與GTF中的基因注釋信息進行比對（Pertea and Pertea，2020），以StringTie的選項“-e-B-P20”估計轉(zhuǎn)錄本豐度，豐度結(jié)果以“.balltown”結(jié)尾的文件呈現(xiàn)，并使用prepDE.py比較這些文件。

1.4 差異表達基因篩選

采用FPKM法計算各樣本中所有基因的表達量，通過DESeq2篩選不同組間的差異表達基因。篩選原則：錯誤發(fā)現(xiàn)率（FDR）＜0.05和|log2FC|＞1的基因即被認為是不同生長階段間的差異表達基因。

1.5 樣本時間序列分析

采用TCseq分析3個生長階段的差異表達基因并進行聚類分析，直觀顯示不同生長階段的不同基因表達模式（Ma et al.，2021）。

1.6 差異表達基因功能富集分析

通過clusterProfiler分別進行差異表達基因的GO功能注釋分析和KEGG信號通路富集分析（Yu et al.，2012），F(xiàn)DR＜0.05即為顯著富集。

2.1 不同生長階段日本醫(yī)蛭的發(fā)育性狀及組織學觀察結(jié)果

經(jīng)4%多聚甲醛固定的日本醫(yī)蛭幼蟲（L）、青年蟲（Y）和成蟲（A），如圖1-A所示。計算不同生長階段日本醫(yī)蛭的體質(zhì)量，其中，幼蟲體質(zhì)量為0.136 g、青年蟲體質(zhì)量為1.184 g、成年蟲體質(zhì)量為1.524 g，不同生長階段日本醫(yī)蛭的體質(zhì)量存在極顯著差異（P＜0.01）（圖1-B）。此外，不同生長階段日本醫(yī)蛭的石蠟切片觀察結(jié)果（圖1-C）顯示，日本醫(yī)蛭成年蟲期切片中的肌纖維明顯多于幼蟲期和青年蟲期。

圖1 日本醫(yī)蛭的表型數(shù)據(jù)及石蠟切片觀察結(jié)果Fig.1 Phenotypic data and paraffin section of H.nipponia

2.2 RNA測序分析結(jié)果

構(gòu)建24個日本醫(yī)蛭cDNA文庫，將不同生長階段不同部位的數(shù)據(jù)進行組合，從日本醫(yī)蛭幼蟲期、青年蟲期和成年蟲期各獲得3組樣品（HNL1~HNL3、HNY1~HNY3和HNA1~HNA3），共產(chǎn)生541124257條原始序列（Raw reads）；經(jīng)Fastp質(zhì)量控制后，獲得517805628條有效序列（Clean reads）。所有樣本的GC含量十分接近，位于41%~45%之間（表2）。以上結(jié)果表明測序數(shù)據(jù)質(zhì)量較好，可滿足后續(xù)分析的要求。

表2 日本醫(yī)蛭轉(zhuǎn)錄組測序數(shù)據(jù)統(tǒng)計結(jié)果Table 2 Summary statistics of H.nipponia transcriptome sequencing

RNA-Seq共組裝出20430個mRNAs轉(zhuǎn)錄本，其中88個為新發(fā)現(xiàn)的mRNAs轉(zhuǎn)錄本。不同樣本的基因表達情況如圖2-A所示。主成分分析（PCA）結(jié)果表明，9個樣品的基因表達量數(shù)據(jù)分布基本一致（圖2-B），說明不同分組樣本間可較清晰地劃分，組內(nèi)相關(guān)性也較高，可為后續(xù)分析提供可靠數(shù)據(jù)。

圖2 mRNAs表達測序分析結(jié)果Fig.2 mRNAs expression sequencing analysis

2.3 差異表達基因及樣本時間序列分析結(jié)果

根據(jù)FDR＜0.05和|log2FC|＞1的差異表達基因篩選標準，對日本醫(yī)蛭不同生長階段進行兩兩對比，結(jié)果共獲得1348個差異表達基因。其中，在幼蟲期—青年蟲期（HNL vs HNY）中獲得238個差異表達基因，在幼蟲期—成蟲期（HNL vs HNA）中獲得976個差異表達基因，在青年蟲期—成蟲期（HNY vs HNA）中獲得763個差異表達基因（圖3）。可見，成蟲期較其他2個生長期的基因表達差異明顯，表現(xiàn)活躍；而幼蟲和青年蟲間的差異相對較小，表現(xiàn)較穩(wěn)定，說明生長發(fā)育階段對日本醫(yī)蛭基因表達影響顯著。進一步分析差異表達基因的表達趨勢，結(jié)果（圖4）發(fā)現(xiàn)基因表達趨勢表現(xiàn)為：3個樣本（Cluster 2、Cluster 4和Cluster 6）持續(xù)下降、2個樣本（Cluster 3和Cluster 5）持續(xù)上升、3個樣本（Cluster 1、Cluster 7和Cluster 8）青年期高表達和1個樣本（Cluster 9）青年期低表達。

圖3 日本醫(yī)蛭差異表達基因統(tǒng)計結(jié)果Fig.3 Statistics of H.nipponia DEGs

圖4 日本醫(yī)蛭差異表達基因樣本時間序列分析結(jié)果Fig.4 DEGs sample time series analysis of H.nipponia

2.4 差異表達基因樣本時間序列GO功能注釋和KEGG信號通路富集分析結(jié)果

為揭示調(diào)控日本醫(yī)蛭發(fā)育的分子機制，對其上調(diào)趨勢的Cluster 4（含329個差異表達基因）和下調(diào)趨勢的Cluster 5（含176個差異表達基因）進行GO功能注釋分析，結(jié)果（圖5）表明，Cluster 4主要注釋在肌球蛋白復(fù)合體、運動活動、肌動蛋白絲結(jié)合、蛋白運輸及金屬酶等功能條目上，Cluster 5主要注釋于離子轉(zhuǎn)運、信號傳導(dǎo)及氧運輸?shù)裙δ軛l目上。KEGG信號通路富集分析結(jié)果（圖6）顯示，差異表達基因樣本時間序列主要富集在氨基酸代謝、苯丙氨酸、酪氨酸和色氨酸生物合成及原核生物碳固定等通路上，但富集結(jié)果顯著性不高（FDR＞0.05）。

圖5 差異表達基因樣本時間序列GO功能注釋分析結(jié)果Fig.5 GO functional annotation analysis of sample time series of DEGs

圖6 差異表達基因樣本時間序列KEGG信號通路富集分析結(jié)果Fig.6 KEGG signal pathway enrichment analysis of sample time series of DEGs

2.5 不同生長階段差異表達基因GO功能注釋和KEGG信號通路富集分析結(jié)果

對不同生長階段差異表達基因進行GO功能注釋分析和KEGG信號通路富集分析，結(jié)果發(fā)現(xiàn)，幼蟲期與成蟲期相比，上調(diào)差異表達基因顯著注釋在肌球蛋白復(fù)合體、細胞骨架、金屬肽酶活性和金屬內(nèi)肽酶活性等GO功能條目上（圖7-A），下調(diào)差異表達基因主要注釋在幾丁質(zhì)代謝過程、對傷害的反應(yīng)、氧結(jié)合、血紅素結(jié)合和離子運輸?shù)菺O功能條目上（圖7-B）；上調(diào)差異表達基因主要富集于氨基酸代謝、糖胺聚糖生物合成—硫酸軟骨素/硫酸皮膚素等KEGG信號通路上（圖7-C），但富集不顯著（FDR＞0.05），下調(diào)差異表達基因則顯著富集在丙烯酰胺轉(zhuǎn)移酶、卟啉與葉綠素代謝、苯丙氨酸/酪氨酸/色氨酸生物合成等KEGG信號通路上（圖7-D）。此外，其他不同生長階段的差異表達基因的GO功能注釋和KEGG信號通路富集分析結(jié)果顯示，上調(diào)的代謝通路主要富集在能量代謝和肌肉發(fā)育方面，而下調(diào)的通路富集在代謝和信號傳導(dǎo)方面（表3）。這可能是日本醫(yī)蛭幼蟲在生長過程中大量攝食并通過消化、代謝和合成物質(zhì)為自身發(fā)育提供能量，而成年日本醫(yī)蛭發(fā)育成熟后其脫皮現(xiàn)象已完成，與其生長發(fā)育相關(guān)的幾丁質(zhì)代謝過程、生物體必需氨基酸合成通路均下調(diào)，尤其是成年日本醫(yī)蛭已適應(yīng)環(huán)境，對于環(huán)境刺激的敏感度也有所降低。

表3 差異表達基因顯著富集的KEGG信號通路Table 3 KEGG signaling pathway with significant enrichment of DEGs

圖7 日本醫(yī)蛭幼蟲期—成年期差異表達基因的功能富集分析結(jié)果Fig.7 The DEGs enrichment analysis in HNL vs HNA periods of H.nipponia

2.6 日本醫(yī)蛭發(fā)育性狀相關(guān)的潛在基因

進一步對差異表達基因的富集途徑分析發(fā)現(xiàn)，GO功能注釋獲得肌球蛋白復(fù)合體、細胞骨架等組分中的差異表達基因在日本醫(yī)蛭的生長發(fā)育過程中（HNL→HNY→HNA）呈上調(diào)趨勢（圖8-A）；而表面蛋白（SPA17）、IQ motif containing G（IQCG）、ACTB、5-羥色胺（5-HT）和HOX基因家族等同源基因可能參與日本醫(yī)蛭的性腺發(fā)育或產(chǎn)卵調(diào)控（圖8-B）；參與神經(jīng)發(fā)育相關(guān)的微管蛋白、細絲蛋白A（FLNA）、原鈣黏蛋白（PCDHGA10）、溶質(zhì)載體家族5成員8（SLC5A8）和NK型同源盒基因轉(zhuǎn)錄因子2.8（NKX2-8）等基因在各生長階段均有表達，且以在幼蟲期的表達量最高（圖8-C）。本研究還發(fā)現(xiàn)，與抗凝相關(guān)的蛋白基因在青年蟲中的表達量較高，包括水蛭素（Hirudin）、抗蛋白酶（ANTA）、金屬蛋白酶家族（ADAMTS）、Therostasin和血小板糖蛋白抑制劑（DECO）等（圖8-D）。

圖8 日本醫(yī)蛭生長發(fā)育相關(guān)差異表達基因的表達情況Fig.8 The expression of DEGs related to growth and development of H.nipponia

3 討論

本研究對不同所長階段的日本醫(yī)蛭進行轉(zhuǎn)錄組學分析，以期為研究日本醫(yī)蛭發(fā)育過程中的生理機制及提高其人工養(yǎng)殖效益提供新的理論依據(jù)。在3個生長階段中共發(fā)現(xiàn)1348個差異表達基因，且差異表達基因數(shù)量隨著日本醫(yī)蛭的發(fā)育而逐漸減少，可能與3月齡青年蟲已發(fā)育至成蟲有關(guān)。差異表達基因的GO功能注釋分析及KEGG信號通路富集分析結(jié)果表明，肌球蛋白復(fù)合體、細胞骨架和細胞生長發(fā)育相關(guān)信號通路呈上調(diào)趨勢，是由于幼蟲和青年蟲時期日本醫(yī)蛭生長速度快，需大量攝食以獲得能量，與非生殖階段（攝食階段）成蟲期相比觀察到更積極攝食行為的結(jié)論（Xu et al.，2021）一致；下調(diào)差異表達基因則主要與信號轉(zhuǎn)導(dǎo)、苯丙氨酸、酪氨酸和色氨酸生物合成等信號通路相關(guān)，究其原因是隨著日本醫(yī)蛭的發(fā)育，其對食物的需求度及代謝水平均逐漸下降。

已有研究表明，動物的比生長速率在不同生長階段存在顯著變化，說明潛在的肌肉生長調(diào)節(jié)機制可能不同（Zhao et al.，2021）。本研究發(fā)現(xiàn)，成年日本醫(yī)蛭的肌纖維數(shù)量明顯多于幼蟲和青年蟲，GO功能注釋分析也顯示HNL vs HNA和HNY vs HNA的上調(diào)差異表達基因主要富集在肌球蛋白復(fù)合體、細胞骨架等條目上。其中，副肌球蛋白是環(huán)節(jié)動物的肌肉組成蛋白，在肌肉緊張收縮機制中發(fā)揮重要作用（Sonobe et al.，2016；Li et al.，2021）。原肌球蛋白是肌肉細肌絲中與肌動蛋白結(jié)合的蛋白，除了完成基本的收縮功能外，還參與完成體內(nèi)多種生命活動，如胞質(zhì)分裂及信號轉(zhuǎn)導(dǎo)等；通過與肌動蛋白的細肌絲螺旋形結(jié)合存在，而調(diào)節(jié)與肌球蛋白間的相互作用（Pavadai et al.，2020）。ACTB是肌動蛋白家族的重要成員之一，幾乎參與所有形式的細胞和細胞器運動，在生長發(fā)育、生物調(diào)控及修復(fù)過程中發(fā)揮重要作用（Hammer et al.，2019）。肌球蛋白由MYL和MYH組成，其中，MYH家族是調(diào)節(jié)肌肉發(fā)育的關(guān)鍵成員之一，主要通過調(diào)節(jié)細胞增殖、分化來影響肌肉的生成（Liu et al.，2020；Hill et al.，2021），這些基因均在成年日本醫(yī)蛭中高表達。

抗凝血能力是水蛭研究的熱點之一。本研究發(fā)現(xiàn)，在日本醫(yī)蛭體內(nèi)除了水蛭素以外，抗凝血基因ANTA和DECO也有表達。水蛭素可抑制胰蛋白酶、糜蛋白酶、組織蛋白酶G和組織激肽釋放酶，但不抑制凝血因子Xa活性（王斌等，2017；Lu et al.，2018）。相反，ANTA是凝血因子Xa的有效抑制劑。DECO能抑制纖維蛋白原表達，與糖蛋白IIb-IIIa復(fù)合物上的血小板受體相互作用，防止血液凝結(jié)（Kwak et al.，2019）?？鼓蛟谌毡踞t(yī)蛭青年期的表達量高于成年，該結(jié)論為青年日本醫(yī)蛭的藥用提供了理論依據(jù)。

水蛭是典型的雌雄同體，雄性腺優(yōu)先于雌性腺發(fā)育成熟。與其他品種水蛭相比，日本醫(yī)蛭的繁殖性能較低（馬燕，2016）。本研究通過轉(zhuǎn)錄組測序分析不同生長發(fā)育階段的差異表達基因，篩選出與性腺發(fā)育或卵發(fā)生相關(guān)的候選基因或蛋白。其中，SP17是一類與精子發(fā)生和形成、精卵識別及頂體反應(yīng)等受精活動密切相關(guān)的蛋白（Minhas et al.，2020）。余米等（2021a）研究證實，在18℃下SP基因在日本水蛭體內(nèi)高表達。IQCG也參與精子的發(fā)生過程，是纖毛/鞭毛運動的關(guān)鍵調(diào)節(jié)器，與精子運動有關(guān)（Harris et al.，2014；Pan et al.，2016）。關(guān)于卵巢發(fā)育的研究表明，5-HT在性腺成熟和連續(xù)產(chǎn)卵方面發(fā)揮多種生物學作用，是通過G蛋白偶聯(lián)受體超家族5-HT受體發(fā)揮作用，調(diào)節(jié)無脊椎動物和脊椎動物的卵母細胞發(fā)育與成熟（Song et al.，2016）。HOX基因是一個調(diào)節(jié)分子家族且編碼高度保守的轉(zhuǎn)錄因子，其表達受性類固醇的調(diào)控，在生殖道發(fā)育、子宮內(nèi)膜周期生長和胚胎植入中發(fā)揮重要作用（Du and Taylor，2015；Kondo et al.，2017）。這些基因在日本醫(yī)蛭成年期顯著上調(diào)，可能與其繁殖性能相關(guān)。

日本醫(yī)蛭是研究神經(jīng)再生機制的理想模型系統(tǒng)，具有在整個生命周期內(nèi)再生中樞神經(jīng)系統(tǒng)的能力。此外，日本醫(yī)蛭中樞神經(jīng)元持續(xù)擴張，即其生長和增加突觸耦合的機制可能永遠不會被完全關(guān)閉（Meriaux et al.，2011；Cohen et al.，2020）。本研究發(fā)現(xiàn)，微管蛋白（刁磊等，2019）、FLNA（周翔宇，2021）及軸突生長誘向因子（Netrin）（孟凡偉等，2003）等與神經(jīng)發(fā)育相關(guān)的基因在3個生長階段均有表達，其中，F(xiàn)LNA、PCDHGA10（Schoch et al.，2017）和NKX2-8在幼蟲期高表達，故可作為篩選日本醫(yī)蛭的候選分子遺傳標記。

4 結(jié)論

在幼蟲期—成蟲期（HNL vs HNA）獲得的976個差異表達基因中，上調(diào)差異表達基因主要與肌肉發(fā)育相關(guān)，而與性腺發(fā)育或產(chǎn)卵調(diào)節(jié)及神經(jīng)發(fā)育相關(guān)的差異表達基因具有明顯周期特異性，其表達變化均與日本醫(yī)蛭的生長發(fā)育有關(guān)。水蛭素和抗凝相關(guān)基因在青年期高表達，因此開展水蛭抗凝血能力研究時宜選用青年日本醫(yī)蛭。