基于轉錄組測序的方斑東風螺單核苷酸多態性位點挖掘及功能注釋

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

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

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

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

隨著高通量測序技術的快速發展，轉錄組測序成本大幅下降，利用比較轉錄組學方法和序列比對識別大量的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位點進行多態特征分析，但未見關于方斑東風螺SNP位點的研究。通過轉錄組測序技術，結合所測物種基因組信息，更易找到與目標性狀相關的SNP[27-28]。本研究通過對方斑東風螺轉錄組的深度測序分析，篩選出大量SNP位點，并對這些SNP所在基因進行功能注釋，為方斑東風螺的抗病及育種研究提供基礎數據。

1 材料和方法

1.1 材料

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

1.2 轉錄組測序數據

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

1.3 SNP位點的檢測

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

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

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

2 結果與分析

2.1 轉錄組測序及序列組裝

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

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

2.2 轉錄組數據及SNP位點數據

利用在對照組轉錄組（BLANK）和實驗組轉錄組（LPS-4h，LPS-8h）獲得的數據，經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位點分析統計發現，BLANK 轉錄組中含1 個SNP位點的unigene 有9 802 條（26.39%），含2～10 個SNP位點的unigene 21 224 條（57.15%），含10 個以上SNP位點的unigene 6 110 條（16.45%）；LPS-4h 轉錄組中含有1 個SNP位點的unigene 有9 721 條（26.22%），含2～10 個SNP位點的unigene 21 224條（57.19%），含10 個以上SNP位點的unigene 6 110條（16.59%）；LPS-8h 轉錄組中含1 個SNP位點的unigene 有9 580 條（26.13%），含2～10 個SNP位點的unigene 20 921 條（57.07%），含10 個以上SNP位點的unigene 6 156 條（16.79%）（圖1）。

表1 測序數據質量分析Table 1 Quality analysis of sequencing data

表2 單基因簇統計分析Table 2 Statistics analysis of Unigenes

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

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

圖1 SNP的分布統計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 結果顯示，共有16 891 條SNP-unigenes匹配相應的COG 注釋信息；根據功能信息可分為26 類，其中“僅通用功能預測”和“信號轉導機制”類最多，分別包含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 子集中，其中以“內吞作用”子集中富集的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 免疫防御相關SNP-unigene的富集及特異SNP位點

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

2.5 差異表達免疫防御相關基因SNP位點

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

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

表5 BLANK、LPS-4h、LPS-8h 轉錄組特異分布免疫防御相關基因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 檢測有助于了解不同個體和群體對外部刺激的反應[30]。因此，開發與方斑東風螺免疫相關的功能基因，挖掘與免疫性狀連鎖的分子標記是全面了解方斑東風螺免疫反應，開展分子標記輔助育種的基礎。本研究基于Hiseq-2000測序技術分別從BLANK、LPS-4h、LPS-8h 轉錄組中獲得224 055、225 287、224 440 個SNP位點。SNP堿基替換類型分為轉換 (A-G、T-C) 和顛換(A-T、C-G、A-C、T-G) 兩類。理論上，顛換與轉換比例應為1∶2[31]。本研究中，轉錄組的顛換與轉換SNP數比值約為1∶1.72，與先前研究中堿基轉換類型發生頻率高于顛換類型類似[31-33]；在6 種核苷酸變異類型中，以A/G、C/T 轉換最多，約占6種核苷酸變異類型的60%，與其他物種的堿基替換頻率類似[34-37]。

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

4 結語

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