余 米，封孝蘭，周 夢，張承露，曹 敏

(重慶市藥物種植研究所，重慶 408435)

日本醫(yī)蛭(Hirudo nipponia)俗稱水蛭、日本醫(yī)水蛭或稻田醫(yī)蛭，屬醫(yī)蛭行亞目(Hirudiniformes Caballero)醫(yī)蛭科(Hirudinidae Whitman)醫(yī)蛭屬(Hirudo)，在中國南北多省均有分布[1]。由于其含有高效抗凝血物質(zhì)-水蛭素，是心腦血管疾病藥物的主要動物源材料。近年來，隨著中藥材市場和藥企需求的增加，過度捕撈加之天然水域水質(zhì)的改變，使得天然水域日本醫(yī)蛭的資源量減少和個體趨于小型化[2]，各地對日本醫(yī)蛭的養(yǎng)殖逐漸興起。

高通量測序技術(shù)具有成本低、不需要預(yù)先得到物種基因信息和獲得的數(shù)據(jù)豐富等優(yōu)點，可用于研究不同蛭類對不同環(huán)境的適應(yīng)內(nèi)在機制[3]、寬體金線蛭在飼喂前后抗凝血因子的基因表達[4]以及日本醫(yī)蛭在饑餓脅迫后唾液腺基因表達水平的差異[5]。以上研究都集中在蛭類的進化和抗凝血活性方面，有關(guān)日本醫(yī)蛭不同組織的基因轉(zhuǎn)錄和SSR 分析未見報道。溫度是影響日本醫(yī)蛭生理的主要因素，18 ℃是日本醫(yī)蛭繁殖的適宜溫度，24 ℃是其生長的適宜溫度[2]，2 種溫度下日本醫(yī)蛭不同組織的轉(zhuǎn)錄組與SSR 信息分析能為篩選日本醫(yī)蛭性別特異分子標(biāo)記、分子標(biāo)記輔助育種以及提高抗凝血活性等提供基礎(chǔ)數(shù)據(jù)。

本試驗以日本醫(yī)蛭肌肉(M)、生殖環(huán)帶(C)和頭部(H) 3 種組織為研究對象，采用新一代高通量測序技術(shù)進行轉(zhuǎn)錄組測序，通過數(shù)據(jù)組裝、組織表達量比較、差異表達基因注釋和SSR 序列分析，以期為日本醫(yī)蛭繁殖和抗凝血活性的分子機理以及分子標(biāo)記輔助育種提供理論依據(jù)。

1 材料與方法

1.1 樣品采集

采用來自重慶市藥物種植研究所水蛭繁育基地的40 尾日本醫(yī)蛭，體質(zhì)量為1.425～2.304 g，體長9.80～10.52 cm。2019 年11 月，開始分別置于24 和18 ℃恒溫培養(yǎng)箱內(nèi)饑餓處理28 d，分別剪取其肌肉(M)、環(huán)帶(C)和頭部(H)組織，浸泡在RNA-hold 保存液中備用，24 ℃處理的日本醫(yī)蛭不同組織簡寫為：H24、C24 和M24；18 ℃處理的日本醫(yī)蛭不同組織簡寫為：H18、C18 和M18。

1.2 總RNA 提取

利用Trizol 法提取不同組織RNA。采用瓊脂糖凝膠電泳分析樣品RNA 完整性及是否存在DNA 污染，利用NanoPhotometer spectrophotometer 檢測RNA 純度(OD260/OD280及OD260/OD230比值)；利用Agilent 2 100 bioanalyzer 精確檢測RNA 完整性。將得到的RNA 溶液置于-80 ℃保存或立刻用于后續(xù)試驗。

1.3 建庫測序和拼接組裝

日本醫(yī)蛭轉(zhuǎn)錄組測序和拼接組裝委托北京諾禾致源科技股份有限公司完成，建庫中使用的建庫試劑盒為Illumina 的 NEBNext?UltraTMRNA Library Prep Kit，建庫完成后使用 Agilent 2 100 bioanalyzer 對文庫的 insert size 進行檢測，insert size 符合預(yù)期后，采用qRT-PCR 對文庫有效濃度進行準(zhǔn)確定量(文庫有效濃度高于 2 nmol/L)。獲得clean reads 后，采用Trinity[6]對clean reads 進行無參拼接。

1.4 功能注釋

為獲得相應(yīng)基因功能的注釋信息，使用BLAST 軟件(2.2.28+版本)[7]在線與非冗余蛋白庫(non-redundant protein sequence database，NR)、核酸序列數(shù)據(jù)庫(nucleotide sequence database，Nt)、Swiss-prot 蛋白序列數(shù)據(jù)(manually annotated and reviewed protein sequence database，Swissprot)和真核生物蛋白質(zhì)同源簇數(shù)據(jù)庫(eukaryotic ortholog groups，KOG)比對，E-value 值均為1e-5；使用KAAS (r140224)軟件[8]與日本京都基因和基因組百科全書(kyoto encyclopedia of genes and genomes，KEGG)比對，E-value 值均為1e-10；使用Hmmscan (HMMER 3)軟件與蛋白質(zhì)家族的集合數(shù)據(jù)庫(protein family，Pfam)[9]比對，E-value值均為0.01；使用Blast2go (b2g4pipe_v2.5)[10]基因本體數(shù)據(jù)庫(gene ontology，GO)比對，E-value值為1.0e-6。

1.5 SSR 標(biāo)記篩選

采用MISA 1.0 版(默認(rèn)參數(shù))對應(yīng)的各個unit size 的最少重復(fù)次數(shù)分別為：1～10、2～6、3～5、4～5、5～5、6～5 (例如1～10，單核苷酸為重復(fù)單位時，其重復(fù)數(shù)至少為10 才可被檢測到；2～6，以雙核苷酸為重復(fù)單位時，其最少重復(fù)數(shù)為6)對unigene 進行SSR 檢測。采用Primer3 2.3.5 版(默認(rèn)參數(shù))進行 SSR 引物設(shè)計。

1.6 差異基因表達分析

將 Trinity 拼接得到的轉(zhuǎn)錄組作為參考序列(RefSeq)，將每個樣品的clean reads 往RefSeq 上做mapping，分別過濾掉低質(zhì)量的reads、帶接頭污染的reads 以及含N 的reads。比對過程中我們采用了RSEM軟件[11]，RSEM 中使用到的bowtie2參數(shù)mismatch 0 (bowtie2 默認(rèn)參數(shù))。

2 結(jié)果與分析

2.1 數(shù)據(jù)質(zhì)量評估及序列組裝分析

經(jīng)過測序質(zhì)量控制，共得到125.08 G 的clean data，過濾后的堿基分布平穩(wěn)，平均GC 含量為42.58%，不同溫度下日本醫(yī)蛭環(huán)帶組織樣品Q30 堿基百分比均不小于92.01% (表1)。

將過濾后的reads 片段進行聚類及拼接組裝(表2)，共得到202 570 條transcript 和52 206 條unigene，transcript 的 N50 和N90 (拼接轉(zhuǎn)錄本按照長度從長到短排序，并進行長度累加，當(dāng)累加值不小于unigene 總長50%/90%的數(shù)值就是N50/N90)分別為2 894 和820，unigene 的N50 和N90分別為2 483 和520，組裝完整性較高。transcript和unigene 的序列長度分布于300～2 000 bp 以上，其中＞2 000 bp 區(qū)間的transcript 分布最多(為65 746，占32.5%)，unigene 在300～500 bp 分布最多(為17 231，占33%) (圖1)。transcript 和 unigene 序列的平均長度分別為1 826 和1 390 bp，最大和最小長度均為30 853 bp 和301 bp。

2.2 功能基因注釋與表達量分布

將unigene 序列與7 大數(shù)據(jù)庫進行比對，得出了unigene 在每個數(shù)據(jù)庫中的功能注釋信息。表3 顯示：有56.22%的unigene 至少被1 個數(shù)據(jù)庫注釋。Pfam 和GO 數(shù)據(jù)庫均注釋到22 984條，占unigene 總數(shù)的44.02%，其余依次為NR 數(shù)據(jù)庫(21 967 條，42.07%)、Swiss-Prot 數(shù)據(jù)庫 (16 852，32.27%)、KOG 數(shù)據(jù)庫 (9 846，18.86%)、KO 數(shù)據(jù)庫(6 408，12.27%)和NT 數(shù)據(jù)庫(6 379，12.21%)。

表1 不同溫度下日本醫(yī)蛭各組織測序數(shù)據(jù)統(tǒng)計Tab.1 Statistical data of RNA-seq for H.nipponia at different temperatures

表2 拼接轉(zhuǎn)錄本長度分布Tab.2 Distribution of number and quality of unigenes and transcripts

圖1 轉(zhuǎn)錄本與獨立基因拼接長度和頻數(shù)分布Fig.1 Distribution of splice length and frequency of transcripts and unigene

表3 日本醫(yī)蛭獨立基因注釋數(shù)量統(tǒng)計Tab.3 The statistical number of unigene that were functional annotated in H.nipponia

2.2.1 注釋基因的NR 分類

通過與NR 進行比對注釋，可以獲取與日本醫(yī)蛭轉(zhuǎn)錄組序列具有相似性的近緣物種信息，結(jié)果(圖2)顯示：日本醫(yī)蛭與淡水澤蛭、海蠕蟲、鴨嘴舌形貝、舞行波豆蟲和鞭蟲具有較高的序列同源性，其中與淡水澤蛭 (H.robusta)具有最高的序列相似性(42%)。

圖2 基于NR 數(shù)據(jù)庫比對的同源基因物種來源分布Fig.2 Species distribution of the top BLASTx hits against the NR database for unigene in H.nipponia

2.2.2 注釋基因的GO 分類

對unigene 序列進行GO 注釋之后，按照GO 3 個基本大類(生物學(xué)過程、細(xì)胞組分和分子功能)將注釋成功的全部基因按照進一步細(xì)分的功能進行聚類(圖3)?；谛蛄型葱浴⒈环殖? 大主要類別共有56 個功能群。在GO 分類的生物學(xué)過程中，注釋基因數(shù)量在二級分類排在前3 位的分別是細(xì)胞過程(13 933 條)、新陳代謝過程(12 017 條)和單有機體過程(10 944 條)。在細(xì)胞組分分類中，注釋數(shù)量前三的依次為細(xì)胞(8 385 條)、細(xì)胞成分(8 385 條)和細(xì)胞膜(6 300條)。在分子功能分類中，注釋最多和次多的分別為連接功能(13 012 條)與催化活性(9 623 條)。

2.2.3 注釋基因的KOG 分類

KOG 數(shù)據(jù)庫歸為26 個組，將通過KOG 成功注釋的unigene 按照分組進行分類(圖4)，在KOG 數(shù)據(jù)庫中共注釋到9 849 條。在這26 個分類中，聚類數(shù)量最多的為信號轉(zhuǎn)導(dǎo)機制(1 553條，15.77%)，其次是功能預(yù)測(1 365 條，13.86%)，第三是翻譯后修飾、蛋白翻轉(zhuǎn)和分子伴侶(1 109 條，11.26%)，聚類最少的是未修飾蛋白(2 條)。

2.2.4 注釋基因的KEGG 分類

以KEGG 代謝途徑數(shù)據(jù)庫為依據(jù)，對unigene 的注釋結(jié)果進行分類(圖5)。日本醫(yī)蛭的unigene 代謝途徑可分為細(xì)胞過程(A)、環(huán)境信息處理(B)、遺傳信息處理(C)、代謝(D)和有機系統(tǒng)(E) 32 個亞類；其中，涉及unigene 數(shù)量最多的5 個代謝途徑分別為B 分支中的信號轉(zhuǎn)導(dǎo)(969 條)、E 分支的內(nèi)分泌系統(tǒng)(575 條)、C 分支的翻譯(553 條)、A 分支中的運輸與分解代謝(493 條)和細(xì)胞群落(388 條)。

2.3 unigene 的SSR 分析

圖3 日本醫(yī)蛭GO 分類Fig.3 GO classification

圖4 日本醫(yī)蛭KOG 分類Fig.4 KOG classification

SSR 被稱為短串聯(lián)重復(fù)序列或微衛(wèi)星DNA標(biāo)記。對日本醫(yī)蛭轉(zhuǎn)錄組的52 206 條unigene 序列進行SSR 分析，共鑒定出29 926 個SSR，總長度72 566 318 bp。這些SSR 分布于16 094 個序列中，其中有7 188 條含有1 個以上SSR，其中在化學(xué)成分形成中的SSR 數(shù)量有5 664 個。圖6、7 表明：三核苷酸SSR 數(shù)量最多，有22 167 個，占6 種核苷酸重復(fù)類型的74.07%，其次是單核苷酸(4 519 個，15.10%)、四核苷酸(1 860 個，6.22%)、二核苷酸(1 174 個，3.92%)、六核苷酸(126 個，0.42%)和五核苷酸(80 個，0.27%)。

2.4 日本醫(yī)蛭SSR 核苷酸基序重復(fù)類型及其頻率特征分析

圖5 日本醫(yī)蛭的KEGG 分類Fig.5 KEGG classification

從日本醫(yī)蛭轉(zhuǎn)錄組SSR 核苷酸基序類型來看，其中29 926 個SSR 位點包含了92 種重復(fù)基元。由表4 所示：單核苷酸至六核苷酸重復(fù)分別有2、4、10、23、24 和29 種；出現(xiàn)最多的基元是ATC/ATG (9 496 個，31.73%)，其次是AAT/ATT (7 075 個，23.64%)和AAC/GTT (3 083 個，10.30%)。二核苷酸重復(fù)基元以AC/GT、AG/CT 和AT/AT 為主，分別占56.81%、21.98%和20.70%。三核苷酸重復(fù)基元以ATC/ATG、AAT/ATT 和AAC/GTT 為主，分別占 42.84%、31.92% 和13.90%。四核苷酸重復(fù)基元以ATCC/ATGG、AAAT/ATTT 和AATG/ATTC 為主，分別占18.75%、12.5%和11.25%。五核苷酸重復(fù)基元以AATAT/ATATT、AACAT/ATGTT 和AAAAT/ATTTT 為主，分別占35.59%、18.17%和14.95%。六核苷酸重復(fù)基元以AATGAT/ATCATT、ACCATC/ATGGTG 和AATGAT/ATCATT 為主，分別占59.52%、4.76%和4.76%。二核苷酸基元重復(fù)數(shù)主要為6次，三～六核苷酸基元重復(fù)數(shù)均為5 次。

2.5 日本醫(yī)蛭不同部位差異基因表達分析

首先對原始的read count 進行標(biāo)準(zhǔn)化(normalization)，主要是對測序深度的校正；然后通過統(tǒng)計學(xué)模型進行假設(shè)檢驗概率(P-value)的計算，最后進行多重假設(shè)檢驗校正(BH)，得到FDR 值，差異基因篩選標(biāo)準(zhǔn)為：|log2(FoldChange)|＞1 和padj＜0.05。篩選結(jié)果如表5 所示：日本醫(yī)蛭頭部在不同溫度下差異表達434 個基因，上調(diào)211個，下調(diào)223 個；環(huán)帶部位在不同溫度下差異表達591 個基因，上調(diào)191 個，下調(diào)400 個；肌肉組織差異表達578 個，上調(diào)362 個，下調(diào)216 個。

圖6 SSR 重復(fù)單元Fig.6 SSR motif unit

圖7 日本醫(yī)蛭重復(fù)基元種類分布Fig.7 Types and distribution of SSRs

表4 日本醫(yī)蛭SSR 不同核苷酸重復(fù)基元的類型及頻率Tab.4 Types and frequencies of different nucleotide repeat motifs of SSR in H.nipponia

表5 差異表達基因數(shù)Tab.5 Number of differentially expressed genes

如圖8 所示：整體來說，水蛭素基因cluster-9 827.360 6 與cluster-9 827.171 34 在H18 中表達水平最高；基因cluster-9 827.226 42 和cluster-9 827.175 32 在H18、C18 和M18 中表達均上調(diào)；cluster-9 827.891 8 在C18 中的表達水平較高。

圖8 日本醫(yī)蛭不同組織相關(guān)基因表達熱圖Fig.8 Heat map of gene expression in different tissues of H.nipponia

3 討論

日本醫(yī)蛭作為一種重要的心腦血管類藥物的原材料，國內(nèi)在日本醫(yī)蛭的抗凝血方面做了大量的基礎(chǔ)性研究，但在分子育種和代謝層面缺乏相關(guān)的基礎(chǔ)數(shù)據(jù)資料。為了獲取日本醫(yī)蛭不同組織的轉(zhuǎn)錄組信息，本試驗在不同溫度下對其頭部、環(huán)帶和肌肉3 種組織進行高通量轉(zhuǎn)錄組測序，對數(shù)據(jù)進行過濾與質(zhì)控后，采用Trinity 軟件組裝、拼接和轉(zhuǎn)錄本功能注釋，得到125.08 G 轉(zhuǎn)錄組clean data 信息，這些數(shù)據(jù)的積累可為開展日本醫(yī)蛭的繁殖和代謝提供參考資料。從轉(zhuǎn)錄組質(zhì)量來看，共組裝得到平均長度為1 390 bp unigene 共計52 206 條，未發(fā)現(xiàn)＜200 bp 的序列。N50 和N90分別為2 483 和520，Q30 均在92.01%以上，其測序質(zhì)量可靠。

本研究得到的SSR 有29 926 個，unigene 共52 206 條，SSR 出現(xiàn)的頻率是57.32%，遠高于蛭類中的棒紋牛蛭(35.6%)、洞穴山蛭(24.9%)和寬體金線蛭(22.3%)[12]，也高于江鱈(10.27%)[13]和牙鲆(7.95%)[14]等魚類。由于選取的是肌肉、環(huán)帶和頭部不同組織的混合樣本，RNA 來源豐富，組裝的轉(zhuǎn)錄組質(zhì)量較好，所以篩選到的SSR 數(shù)量相對較多，這些大量的序列為發(fā)現(xiàn)日本醫(yī)蛭新基因提供了充足的轉(zhuǎn)錄組序列基源。同時也證實了轉(zhuǎn)錄組測序是開發(fā)SSR 分子標(biāo)記一種有效方法。

本研究所有鑒定到的微衛(wèi)星中，SSR 的核苷酸基元重復(fù)數(shù)都在 5 次以上，主要集中在 5～10次之間(86.27%)，以三核苷酸和單核苷酸重復(fù)SSR 為主，所占比例分別為74.07%和15.10%，六核苷酸和五核苷酸最少，占比分別為0.42%和0.27%。劉洪濤等[15]認(rèn)為富含 A/T 的序列退火溫度較低，有利于 DNA 解鏈以增加其出現(xiàn)的概率，同時由于蛋白質(zhì)密碼子為三堿基造成處于編碼區(qū)的 SSR 三堿基重復(fù)占比最高。另外，本研究結(jié)果顯示：在三核苷酸重復(fù)基元中，ATC 所占比例最高，為9 496 次；單核苷酸重復(fù)基元中，以A 為主，為4 255 次；四核苷酸重復(fù)基元中以ATCC 為主，達到了662 次。這一結(jié)果與王斌等[12]對蛭類的報道存在一定差異，可能是由于物種不同，其SSR 核苷酸基序重復(fù)類型本身就存在差異，例如江鱈以一、二核苷酸基元為主[13]，牙鲆和雙須骨舌魚以二、三核苷酸為主[14,16]。

日本醫(yī)蛭其咽周腺體(位于頭部)分泌的水蛭素，是高效的抗凝血活性物質(zhì)，本次試驗結(jié)果表明：在18 ℃下，日本醫(yī)蛭頭部的水蛭素基因表達明顯(圖 8)，故可在養(yǎng)殖時進行適當(dāng)降溫處理增加日本醫(yī)蛭的水蛭素總體含量。我們在解剖取樣時發(fā)現(xiàn)：經(jīng)過28 d 處理后日本醫(yī)蛭體內(nèi)依然儲藏著大量的未凝固血液，18 ℃相較于24 ℃處理組其儲藏量更大，我們推測降溫后，日本醫(yī)蛭為了防止體內(nèi)血液的凝固而大量分泌水蛭素，而對照組環(huán)境溫度較高，其血液的代謝速度較快，水蛭素的分泌減少。

基因cluster-9 827.226 42 [Sex peptide (SP) family]是性肽家族一員，在18 ℃ 3 種組織中表達均升高，性肽家族在日本醫(yī)蛭體內(nèi)參與雌性交配后的調(diào)節(jié)和信號轉(zhuǎn)導(dǎo)等生物學(xué)過程，果蠅性肽(SP)被認(rèn)為是一種精液成分，可誘導(dǎo)受精雌性的產(chǎn)后反應(yīng)，如抑制再生和刺激產(chǎn)卵，SP 被認(rèn)為在果蠅的性選擇和性沖突中起著中心作用[17-18]。基因cluster-9 827.175 32 (FBXO 家族)在泛素介導(dǎo)的蛋白質(zhì)水解過程中具有底物識別的特性和功能，在本試驗中FBXO8 主要參與ARF 蛋白信號轉(zhuǎn)導(dǎo)、GTPase 酶活性調(diào)節(jié)和DNA 雙鏈斷裂形成等生物學(xué)過程。此外，日本醫(yī)蛭是雌雄同體、異體受精的生物，其在繁殖季節(jié)關(guān)于日本醫(yī)蛭性別的扮演一直未見報道，Dmrt (double sex and mab-3 related transcription factor)是重要的性別決定基因，本研究中發(fā)現(xiàn)日本醫(yī)蛭在18 ℃下性別決定相關(guān)基因cluster-9 827.891 8 在C18 中的表達水平較高，表明18 ℃下日本醫(yī)蛭性別可能發(fā)生轉(zhuǎn)換。由于我們?nèi)迎h(huán)帶部位包括了日本醫(yī)蛭所有生殖腺，故不能確定這些性別決定基因在雌雄生殖腺中的表達情況。但有文獻證實性別決定的假定基因位于性別決定的遺傳圖譜上，具有豐富的小衛(wèi)星標(biāo)記[19-21]。因此，本研究為進一步研究日本醫(yī)蛭繁殖的分子機理和篩選性別特異分子標(biāo)記積累了基礎(chǔ)數(shù)據(jù)，同時有利于抗凝血相關(guān)的分子標(biāo)記的開發(fā)。

4 結(jié)論

不同溫度處理后日本醫(yī)蛭的SSR 發(fā)生頻率較高，抗凝血與繁殖相關(guān)基因(cluster-9 827.360 6、cluster-9 827.226 42 和cluster-9 827.413 等)表達明顯，其頭部和環(huán)帶組織可作為日本醫(yī)蛭篩選抗凝血活性和分子遺傳標(biāo)記的供試組織。