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

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

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

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

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

1 材料和方法

1.1 材料

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

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

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

1.3 SNP位點(diǎn)的檢測(cè)

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

1.4 SNP位點(diǎn)所在 unigene 的注釋與功能分析

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

2 結(jié)果與分析

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

用Illumina 高通量測(cè)序平臺(tái)對(duì)方斑東風(fēng)螺足組織進(jìn)行轉(zhuǎn)錄組測(cè)序，對(duì)原始數(shù)據(jù)進(jìn)行嚴(yán)格的質(zhì)量控制，經(jīng)過(guò)濾后在對(duì)照組轉(zhuǎn)錄組（BLANK）和實(shí)驗(yàn)組轉(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）。說(shuō)明方斑東風(fēng)螺測(cè)序質(zhì)量較高，轉(zhuǎn)錄組數(shù)據(jù)可用于后續(xù)分析。

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

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

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

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

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

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

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

圖1 SNP的分布統(tǒng)計(jì)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ù)測(cè)”和“信號(hào)轉(zhuǎn)導(dǎo)機(jī)制”類最多，分別包含2 848 和2 814 條SNP-unigenes（圖3A）。GO 分析表明，共有4 682 條SNP-unigenes 匹配到GO 條目(GO term)，GO 條目包含生物過(guò)程、細(xì)胞組分及分子功能的42 個(gè)亞類，分別在代謝過(guò)程、催化活性和細(xì)胞部分類中最為富集 (圖3B)。KEGG富集分析顯示，共有5 866 條SNP-unigenes 富集到298 個(gè)KEGG 子集中，其中以“內(nèi)吞作用”子集中富集的SNP-Unigene 最多，共計(jì)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 信號(hào)通路富集分析Fig.4 Kyoto encyclopedia of genes and genomes (KEGG) classification of SNP-unigenes

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

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

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

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

3 討論

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

4 結(jié)語(yǔ)

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