王 菁，劉付柏，許尤厚，黃寶松，王忠良

（1.廣東海洋大學(xué)水產(chǎn)學(xué)院，廣東 湛江 524088；2.廣西北部灣海洋生物多樣性養(yǎng)護重點實驗室，廣西 欽州 535000）

方斑東風(fēng)螺 (Babylonia areolata) 俗稱花螺，隸屬于軟體動物門腹足綱蛾螺目，有生長速度快、養(yǎng)殖周期短、肉味鮮美、軟體部不飽和脂肪酸含量豐富、經(jīng)濟價值高、便于運輸?shù)葍?yōu)點[1-3]。隨著方斑東風(fēng)螺養(yǎng)殖的快速發(fā)展，養(yǎng)殖過程中生長緩慢、病害暴發(fā)等問題日益突出，須通過遺傳改良、病害防控、飼料營養(yǎng)優(yōu)化等解決，其中遺傳改良對促進東風(fēng)螺養(yǎng)殖產(chǎn)業(yè)健康持續(xù)發(fā)展有現(xiàn)實意義[4-5]。

單核苷酸多態(tài)性(Single nucleotide polymorphism,SNP)是指因單個核苷酸變異引起的DNA 序列多態(tài)性，是一種常見的基因突變[6-7]。SNP標記為第3 代分子標記技術(shù)，與限制性片段長度多態(tài)（restriction fragment length polymerphisms，RFLP）、微衛(wèi)星多態(tài)(microsatellite polymorphisms)相比，在基因組中位點豐富，代表性強，遺傳穩(wěn)定，可自動檢測，有低成本、高效率的優(yōu)點[8-10]，因而廣泛應(yīng)用于水產(chǎn)動物研究[11-18]。

隨著高通量測序技術(shù)的快速發(fā)展，轉(zhuǎn)錄組測序成本大幅下降，利用比較轉(zhuǎn)錄組學(xué)方法和序列比對識別大量的SNP位點已逐漸成為一種趨勢。目前，已對櫛孔扇貝(Chlamys farreri)[19]、馬氏珠母貝(Pinctada fucata)[20]、波紋唇魚(Cheilinus undulatus)[21]、曼氏無針烏賊(Sepiella japonica)[22]、棘頭梅童魚(Collichthys lucidus)[23]、大菱鲆(Scophthalmus maximus)[24]、大口黑鱸(Micropterus salmoides)[25]、凡納濱對蝦(Litopenaeus vannamei)[10,26]等SNP位點進行多態(tài)特征分析，但未見關(guān)于方斑東風(fēng)螺SNP位點的研究。通過轉(zhuǎn)錄組測序技術(shù)，結(jié)合所測物種基因組信息，更易找到與目標性狀相關(guān)的SNP[27-28]。本研究通過對方斑東風(fēng)螺轉(zhuǎn)錄組的深度測序分析，篩選出大量SNP位點，并對這些SNP所在基因進行功能注釋，為方斑東風(fēng)螺的抗病及育種研究提供基礎(chǔ)數(shù)據(jù)。

1 材料和方法

1.1 材料

方斑東風(fēng)螺購自廣東省湛江某東風(fēng)螺養(yǎng)殖場，平均體質(zhì)量為25 g，于實驗室海水桶中暫養(yǎng)1 周（80 L，25℃）。實驗時，將方斑東風(fēng)螺分為脂多糖（LPS）注射組（LPS-4h、LPS-8h）和空白對照組。實驗組對方斑東風(fēng)螺閉殼肌注射100 μL 0.5 mg/mL 的LPS懸浮液，對照組注射同體積的磷酸鹽緩沖液（PBS），分別于刺激后4 h、8 h（分別記為LPS-4h、LPS-8h組）采集各組足組織樣品，于液氮中速凍，置-80 ℃下保存?zhèn)溆谩?/p>

1.2 轉(zhuǎn)錄組測序數(shù)據(jù)

樣品委托廣州基迪奧生物科技有限公司使用TRIzol 試劑（Invitrogen，美國）提取總RNA，經(jīng)純度（NanoDrop 2000）和濃度（Agilent 2100）檢測后，同處理組10 個個體樣本合并為1 例樣品進行建庫測序，經(jīng)過Illumina/ Hiseq-2000 高通量轉(zhuǎn)錄組測序，刪除大量低質(zhì)量原始數(shù)據(jù)和適配子等，使用組裝程序Trinity 對轉(zhuǎn)錄組數(shù)據(jù)進行序列組裝(Raw Reads 數(shù)據(jù)的SRA 登錄號為SRP216586)。

1.3 SNP位點的檢測

通過 Call snp 軟件 bcftools (https://github.com/samtools/bcftools) 獲得在不同處理組間有表達差異的SNPs 位點，使用SOAPsnp 對獲得的SNP進行統(tǒng)計和分析。

1.4 SNP位點所在 unigene 的注釋與功能分析

基于所得 SNP-unigenes 序列與 Nr、Nt 及Swiss-prot 數(shù)據(jù)庫比對后的蛋白功能注釋信息，對測序數(shù)據(jù)的KEGG 通路、差異表達基因和SNP進行分析。從轉(zhuǎn)錄組中篩選免疫相關(guān)的KEGG 通路，根據(jù)通路中的免疫基因篩選出差異表達的免疫防御相關(guān)基因，并進行SNP位點分析。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組測序及序列組裝

用Illumina 高通量測序平臺對方斑東風(fēng)螺足組織進行轉(zhuǎn)錄組測序，對原始數(shù)據(jù)進行嚴格的質(zhì)量控制，經(jīng)過濾后在對照組轉(zhuǎn)錄組（BLANK）和實驗組轉(zhuǎn)錄組（LPS-4h、LPS-8h 組）分別獲得56 424 638、50 596 990、49 924 362 條純凈序列，其中GC 比例分別為47.8%、46.15%、46.88%，堿基質(zhì)量值Q30分別為95.40%、95.61%、95.40%（表1）。說明方斑東風(fēng)螺測序質(zhì)量較高，轉(zhuǎn)錄組數(shù)據(jù)可用于后續(xù)分析。

對測序數(shù)據(jù)進行序列組裝，共獲得81 773 條Unigene，總長度為55 763 627 bp，平均長度為681 bp，N50 的長度為1 035 bp (表2)。轉(zhuǎn)錄本和Unigene 的N50 長度均遠大于其平均長度，證明組裝效果較佳。

2.2 轉(zhuǎn)錄組數(shù)據(jù)及SNP位點數(shù)據(jù)

利用在對照組轉(zhuǎn)錄組（BLANK）和實驗組轉(zhuǎn)錄組（LPS-4h，LPS-8h）獲得的數(shù)據(jù)，經(jīng)SOAPsnp軟件檢測，從37 136、37 076、36 657 條unigenes中分別獲得224 055、225 287、224 440 個SNP位點，所有SNP位點中，純合SNP位點107 407 個（BLANK組35 635個，LPS-4h組35 543個，LPS-8h組36 229 個）（表3）。對于Unigene 上的SNP位點分析統(tǒng)計發(fā)現(xiàn)，BLANK 轉(zhuǎn)錄組中含1 個SNP位點的unigene 有9 802 條（26.39%），含2～10 個SNP位點的unigene 21 224 條（57.15%），含10 個以上SNP位點的unigene 6 110 條（16.45%）；LPS-4h 轉(zhuǎn)錄組中含有1 個SNP位點的unigene 有9 721 條（26.22%），含2～10 個SNP位點的unigene 21 224條（57.19%），含10 個以上SNP位點的unigene 6 110條（16.59%）；LPS-8h 轉(zhuǎn)錄組中含1 個SNP位點的unigene 有9 580 條（26.13%），含2～10 個SNP位點的unigene 20 921 條（57.07%），含10 個以上SNP位點的unigene 6 156 條（16.79%）（圖1）。

表1 測序數(shù)據(jù)質(zhì)量分析Table 1 Quality analysis of sequencing data

表2 單基因簇統(tǒng)計分析Table 2 Statistics analysis of Unigenes

表3 SNP位點數(shù)量概況Table 3 SNPidentified in BLANK,LPS-4h and LPS-8h transcriptomes

圖2 表明，BLANK 轉(zhuǎn)錄組的純合SNP位點中，顛換位點 82 447 個(36.80%)，轉(zhuǎn)換位點141 608 個(63.20%)；LPS-4h 轉(zhuǎn)錄組的純合SNP位點中，顛換位點82 878 個 (37.12%)，轉(zhuǎn)換位點142 409 個 (63.78%)；LPS-8h 轉(zhuǎn)錄組的純合SNP位點中，顛換位點82 444 個 (36.73%)，轉(zhuǎn)換位點141 996 個 (63.27%)。在6 種核苷酸的變異類型中，以A/G 轉(zhuǎn)換最多，分別占純合SNP總數(shù)的31.83%（BLANK 組71 316 個）、31.85%（LPS-4h組71 706 個）和31.91%（LPS-8h 組71 630 個）。

圖1 SNP的分布統(tǒng)計Fig.1 SNPdistribution in BLANK,LPS-4h and LPS-8h transcriptomes

圖2 SNP類型分析Fig.2 SNPnumbers of different mutation types in BLANK,LPS-4h and LPS-8h transcriptomes

2.3 SNP-unigene 的注釋與功能

COG 結(jié)果顯示，共有16 891 條SNP-unigenes匹配相應(yīng)的COG 注釋信息；根據(jù)功能信息可分為26 類，其中“僅通用功能預(yù)測”和“信號轉(zhuǎn)導(dǎo)機制”類最多，分別包含2 848 和2 814 條SNP-unigenes（圖3A）。GO 分析表明，共有4 682 條SNP-unigenes 匹配到GO 條目(GO term)，GO 條目包含生物過程、細胞組分及分子功能的42 個亞類，分別在代謝過程、催化活性和細胞部分類中最為富集 (圖3B)。KEGG富集分析顯示，共有5 866 條SNP-unigenes 富集到298 個KEGG 子集中，其中以“內(nèi)吞作用”子集中富集的SNP-Unigene 最多，共計182 條（圖4）。

圖3 SNP-unigenes COG 與GO 功能注釋分析Fig.3 Cluster of orthologous groups (COG) and gene ontology (GO) classification of SNP-unigenes

圖4 SNP-unigenes 的KEGG 信號通路富集分析Fig.4 Kyoto encyclopedia of genes and genomes (KEGG) classification of SNP-unigenes

2.4 免疫防御相關(guān)SNP-unigene的富集及特異SNP位點

根據(jù)KEGG 信號通路的富集分析，篩選到515個 免疫防御相關(guān) SNP-unigenes，注釋到“Autophagy-animal”等19 條與免疫功能相關(guān)的信號通路中，其中以“Autophagy-animal”信號通路中注釋的unigenes 最多(92條)，其次分別為“mTOR signaling pathway”“Wnt signaling pathway”“FoxO signaling pathway”等信號通路(表4)。

2.5 差異表達免疫防御相關(guān)基因SNP位點

根據(jù)轉(zhuǎn)錄組中unigene 的SNP位點分布情況及KEGG 功能注釋信息篩選BLANK、LPS-4h、LPS-8h轉(zhuǎn)錄組中特異分布的免疫防御相關(guān)基因SNP位點，發(fā)現(xiàn)大量的SNP位點存在于免疫防御相關(guān)基因中（表5）?；赗PKM 標準化分析unigene 的表達水平，篩選閾值為q＜0.05，且 |log2（差異倍數(shù)）|＞1，獲得大量差異表達基因。統(tǒng)計這些差異表達基因上的SNP位點，并對其中涉及免疫防御的基因進行KEGG 通路分析。注射LPS 4 h 后，差異表達的免疫防御相關(guān)基因主要參與胞吞作用、IL-17 信號通路、mTOR 信號通路、Wnt 信號通路。注射LPS 8 h 后，存在大量SNP位點的免疫相關(guān)基因主要參與吞噬體、溶酶體通路。比較注射4 h 與8 h，免疫相關(guān)基因則主要參與基座切除修復(fù)、谷胱甘肽代謝、藥物代謝（細胞色素P450）、細胞色素代謝通路（表5）。

表4 免疫防御相關(guān)SNP-unigene 的KEGG 富集分析Table 4 KEGG enrichment analysis of immune-related SNP-unigenes

表5 BLANK、LPS-4h、LPS-8h 轉(zhuǎn)錄組特異分布免疫防御相關(guān)基因SNP位點分析(僅展示前10 行)Table 5 Specific immune-related SNPs identified in BLANK,LPS-4h,LPS-8h transcriptomes (Show only the first 10 rows)

3 討論

在魚類和貝類中，機體對病原體的抵抗力部分受基因控制[29]。開展SNPs 檢測有助于了解不同個體和群體對外部刺激的反應(yīng)[30]。因此，開發(fā)與方斑東風(fēng)螺免疫相關(guān)的功能基因，挖掘與免疫性狀連鎖的分子標記是全面了解方斑東風(fēng)螺免疫反應(yīng)，開展分子標記輔助育種的基礎(chǔ)。本研究基于Hiseq-2000測序技術(shù)分別從BLANK、LPS-4h、LPS-8h 轉(zhuǎn)錄組中獲得224 055、225 287、224 440 個SNP位點。SNP堿基替換類型分為轉(zhuǎn)換 (A-G、T-C) 和顛換(A-T、C-G、A-C、T-G) 兩類。理論上，顛換與轉(zhuǎn)換比例應(yīng)為1∶2[31]。本研究中，轉(zhuǎn)錄組的顛換與轉(zhuǎn)換SNP數(shù)比值約為1∶1.72，與先前研究中堿基轉(zhuǎn)換類型發(fā)生頻率高于顛換類型類似[31-33]；在6 種核苷酸變異類型中，以A/G、C/T 轉(zhuǎn)換最多，約占6種核苷酸變異類型的60%，與其他物種的堿基替換頻率類似[34-37]。

根據(jù)SNPs 所處位置，可將SNPs 分為基因編碼區(qū)SNPs 和基因非編碼區(qū)SNPs[38]，而位于基因調(diào)節(jié)區(qū)的SNP則稱為調(diào)節(jié)SNPs，如其與基因相互作用影響到轉(zhuǎn)錄表達水平，則間接影響RNA 或蛋白質(zhì)的表達數(shù)量，從而引起不同的機體反應(yīng)[39]。本研究中，共有5 866 條SNP-unigenes 富集到298 個KEGG 子集中，以“內(nèi)吞作用”子集中富集的SNP-Unigene 最多。進一步篩選到 515 個SNP-unigenes 注釋到等19 條與免疫功能相關(guān)的信號通路，大部分SNP-unigenes 參與“自噬（動物）”“mTOR 信號通路”“Wnt 信號通路”“FoxO 信號通路”“細胞凋亡”等免疫通路。本研究表明，LPS注射4 h 后，大量差異表達的SNP-unigenes 主要參與胞吞作用、IL-17 信號通路、mTOR 信號通路/ Wnt信號通路。其中，參與胞吞作用的FGFR3基因上存在60 個SNP位點，F(xiàn)GFR 是酪氨酸激酶家族一員，它的激活可能導(dǎo)致腫瘤細胞生長和存活增加[40-41]。Wnt 信號通路是一個高度保守的系統(tǒng)，調(diào)控所有后生動物的復(fù)雜生物學(xué)過程。本研究發(fā)現(xiàn)，Wnt10基因中共存在26 個SNP位點。Wnt10可減少腫瘤細胞間互相黏連，促進腫瘤細胞間轉(zhuǎn)移，進而調(diào)節(jié)干細胞行為、組織穩(wěn)態(tài)和損傷修復(fù)[42]。LPS注射8 h 后，有大量SNP位點的免疫相關(guān)基因主要參與吞噬體、溶酶體通路。參與吞噬體通路的MPO上存在49 個SNP位點，參與溶酶體通路的LITAF基因上存在41 個位點。MPO 是一種血色素蛋白，是活化的中性粒細胞分泌的過氧化物酶類，可能通過各種炎癥反應(yīng)加速動脈粥樣斑塊的氧化，增強巨噬細胞吸收和泡沫細胞形成[43-44]。脂多糖誘導(dǎo)的腫瘤壞死因子 (Lipopolysac-charide-induced TNF-alpha factor，LITAF) 可調(diào)控TNF-α、IL-2 等細胞因子，在無脊椎動物的先天性免疫系統(tǒng)中扮演著重要的介質(zhì)作用[45-46]。本研究LPS-4h 轉(zhuǎn)錄組與LPS-8h 組比較，差異表達的APEX1基因(33 個SNP位點) 和GST基因（32 個SNP位點）主要參與基座切除修復(fù)、谷胱甘肽代謝、藥物代謝 (細胞色素P450)、細胞色素代謝通路。APEX1既有DNA 修復(fù)酶活性，又有氧化還原功能，在多種惡性腫瘤中參與腫瘤的侵襲與轉(zhuǎn)移[47-48]。谷胱甘肽S-轉(zhuǎn)移酶(GST) 作為抗氧化系統(tǒng)中第2 道防線主要成員，催化谷胱甘肽與親電子外源化學(xué)物結(jié)合，將毒性物質(zhì)轉(zhuǎn)化成易排泄的水溶性形式[49]。推測方斑東風(fēng)螺受刺激后，這些SNP-unigenes 可能廣泛參與免疫反應(yīng)[50-52]。

4 結(jié)語

本研究應(yīng)用高通量轉(zhuǎn)錄組測序數(shù)據(jù)實現(xiàn)方斑東風(fēng)螺SNP標記的高效率、大規(guī)模開發(fā)，對存在大量SNP位點的基因進行KEGG 通路分析有助于系統(tǒng)了解方斑東風(fēng)螺在LPS 刺激后免疫反應(yīng)的分子機制，對抗病品系輔助育種研究有重要意義。